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SUMMARY 


The  goal  of  this  project  was  to  improve  and  validate  modeling  tools  for  nonequilibrium 
gas  flows.  Such  flows  are  key  in  many  practical  problems  such  as  shock  layer  radiation  and 
surface  heating  of  atmospheric  re-entry  vehicles  and  spacecraft,  rocket  plume  core  radiation  and 
plume/freestream  interaction,  electric  propulsion  thrusters,  and  low-pressure  plasma-processing 
reactor  design.  Predicting  the  important  physics  of  these  gas  flows  requires  the  ability  to  model 
individual  molecular  collisions  and  simulate  the  energy  transfer  and  chemical  processes  in  a 
computationally  efficient  manner.  The  most  widely  used  technique  for  simulating  gas  flows 
tihrough  computing  individual  collisions  is  known  as  Direct  Simulation  Monte  Carlo  (DSMC). 
Thisproject  has  produced  several  DSMC  codes  for  in-house  simulations  of  immediate  interest. 
In  addition,  this  project  has  undertaken  several  studies  that  address  the  sensitivity  of  DSMC 
results  to  the  various  simplified  energy  transfer  and  chemistry  models  that  are  commonly 
utilized. 

The  first  part  of  this  report  presents  the  results  of  an  investigation  into  the  use  of  shock 
tube  measurements  of  reaction  rate  constants  and  induction  times  to  validate  models  of 
vibrational  relaxation  and  dissociation.  Several  DSMC  vibrational  energy  transfer  models  are 
examined.  It  is  concluded  that  the  strong  coupling  of  vibrational  transfer  and  dissociation  in  a 
shock  tube  precludes  using  the  results  to  validate  the  model  for  one  of  these  processes  without 
having  data  on  the  other.  In  addition,  an  investigation  of  vibrational  favoring  in  air  exchange 
reactions  is  performed,  and  it  is  shown  that  room  temperature  vibrational  disposal  measurements 
of  the  reverse  (exothermic)  reactions  provide  sufficient  information  to  conclude  that  vibrational 
favoring  is  not  important  in  these  cases. 

The  second  part  of  this  report  looks  in  more  detail  at  various  DSMC  models  that  can  be 
used  to  simulate  molecular  dissociation.  It  is  shown  that  reproducing  a  temperature-dependent 
rate  constant  (Arrhenius-type)  is  an  insufficient  validation  for  the  model  to  reasonably  reproduce 
far-from-equilibrium  reacting  flows.  In  the  absence  of  state-specific  data,  quasi-classical 
trajectory  (QCT)  results  for  H2  dissociation  are  used  to  compare  the  qualitative  behavior  of 
different  models  in  terms  of  vibrational  favoring. 

In  the  third  part  of  the  report,  direct  simulation  Monte  Carlo  simulation  results  are 
compared  to  rotational  temperature  measurements  of  nitrogen  in  finite  background  jet  expansion 
flows.  The  influence  of  the  rotational  relaxation  model,  the  influx  boundary  condition,  and  the 
outer  boundary  treatment  on  the  simulation  are  shown  to  have  a  significant  effect  on  the 
predicted  rotational  temperate  shock  structure. 

A  number  of  details  necessary  for  the  proper  DSMC  implementation  of  collision- 
selection  models,  rotational  transfer,  vibrational  transfer,  and  chemistry  have  been  worked  out 
and  are  given  in  the  appendix. 
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PART  I:  INFLUENCE  OF  VIBRATIONAL  NONEQUILIBRIUM  ON  CHEMICALLY 
REACTING  RAREFIED  FLOWS:  TOWARD  EXPERIMENTAL  VERIFICATION  OF 
DSMC  MODELS 


INTRODUCTION 

An  accurate  model  of  chemical  reactions  is  a  key  factor  in  obtaining  useful  predictions  of 
radiation  and  surface  heating  around  re-entry  vehicles  or  spacecraft  and  in  understanding  plasma 
processing  environments,  etc.  These  types  of  flows,  however,  typically  display  significant 
thermal  and  chemical  nonequilibrium.  For  instance,  the  vibrational  temperature  of  the  molecules 
may  often  be  substantially  different  than  the  translational  temperature  of  the  gas.  Indeed, 
molecules  in  these  flows  often  have  a  non-Boltzmann  vibrational  state  distribution  such  that 
there  is  no  true  temperature  associated  with  the  vibrational  degrees  of  freedom.  At  the  same  time, 
the  degree  of  population  of  high  vibrational  levels  is  often  the  controlling  factor  in  the 
probability  of  dissociation  and  endoergic  exchange  reactions.  Recent  work  by  Boyd,  et  al. 1  uses 
Direct  Simulation  Monte  Carlo  (DSMC)  to  predict  density  and  temperature  of  NO  produced  in  a 
re-entry  bow  shock  and  shows  that  different  nonequilibrium  chemistry  models  can  change  the 
predicted  number  density  of  NO,  and  thus  the  predicted  ultraviolet  radiation,  dramatically.  Thus, 
it  is  clear  that  an  accurate  model  for  chemical  reactions  under  conditions  of  vibrational 
nonequilibrium  is  needed,  and  extensive  research  in  this  area  has  taken  place  over  the  past 
several  decades. 

DSMC  calculations,  which  simulate  individual  gas  particle  collisions,  are  well  suited  to 
the  inclusion  of  physically  accurate  chemistry  models  that  are  based  on  the  translational  and 
internal  energy  states  of  each  molecule.  Nonequilibrium  chemistry  models  have  been  the  topic  of 
considerable  research.  Various  models  and  approximations  are  being  implemented  in  DSMC 
codes  for  use  in  predicting  rarefied  flowflelds.  One  could  be  certain  that  a  DSMC  chemistry 
model  is  accurate  for  any  flowfield  if  one  had  a  complete  set  of  state-specific  reaction  and  energy 
transfer  cross  sections  as  a  function  of  collision  parameters  (velocity,  impact  parameter,  etc.)  for 
all  species.  Future  high-level  calculations  and  laboratory  measurements  may  ultimately  provide 
these,  but  in  the  meantime,  cruder  means  of  model  validation  are  sought.  Experimental  data 
against  which  to  compare  nonequilibrium  chemistry  models,  however,  are  quite  difficult  to 
obtain.  Large  shock  tube  facilities  are  one  way  to  try  to  create  high-temperature  nonequilibrium 
gas  flowfields,  but  detailed  molecular  data  from  these  flows  are  rare. 

Preliminary  consideration  by  this  group  of  possible  experiments  that  couple  rarefied 
flowfields  and  chemical  reactions  revealed  the  necessity  for  a  series  of  parametric  calculational 
studies  to  clarify  what  data  will  be  useful.  Many  measurements  of  reacting  flowfields  or  state- 
specific  chemical  reactions  would  be  possible  that  would  not  increase  the  understanding  of  how 
to  properly  implement  a  general  DSMC  chemistry  code  or  distinguish  between  some  of  the 
candidate  models.  For  instance,  Treanor,  et  al.2  have  shown  that  unless  the  degree  of  thermal 
nonequilibrium  in  an  air  shock  is  great  (vibrational  temperature  is  less  than  one-tenth  the 
translational  temperature),  then  differences  between  the  predictions  of  different  chemistry  models 
on  the  production  of  NO  (which  is  the  measured  quantity  in  some  shock  tube  experiments) 
cannot  be  distinguished. 

Shocked  gas  flows  are  typical  of  many  nonequilibrium  flows  in  that  the  translational  and 
rotational  degrees  of  freedom  very  quickly  reach  some  high  temperature,  while  the  vibrational 
degrees  of  freedom  lag  behind  with  a  much  lower  energy.  Because  this  is  the  key  feature  of  die 
nonequilibrium  chemistry  of  many  rarefied  flows,  it  is  possible  to  perform  useful  parametric 
calculational  studies  for  a  highly  simplified  case.  The  temporal  evolution  of  initially  vibrationally 
cold  molecules  in  an  isothermal  bath  (described  in  more  detail  here)  that  is  in  translation-rotation 
equilibrium  displays  many  of  the  same  physical  mechanisms  that  control  the  nonequilibrium 
chemistry  of  shocked  gas  flows. 
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ATOM  EXCHANGE  REACTIONS 

A  wide  variety  of  exoergic  atom  exchange  reactions  are  known  to  produce  products  that 
are  highly  vibrationally  (or  rotationally)  excited,  or  "hot,"  compared  with  a  thermal  or  statistical 
distribution  of  the  available  reaction  energy  into  all  degrees  of  freedom.  Levine,  et  alM  have 
shown  that  the  principle  of  detailed  balance  requires  those  reactions  that  produce  enhanced  post¬ 
reaction  vibrational  excitation  in  the  exoergic  direction  must  also  have  their  reaction  probability 
enhanced  by  vibrational  excitation  of  reactants  in  the  reverse  (endoergic)  direction.  This 
characteristic  will  be  referred  to  as  "vibrational  favoring."  Note  that  any  endoergic  reaction  will 
be  enhanced  by  vibrational  excitation,  merely  because  the  amount  of  total  energy  available  in  the 
^increased  (assuming  that  all  types  of  energy  contribute  to  the  activation  energy). 

Vibrational  favoring"  refers  to  an  increase  in  reaction  probability  beyond  what  would  be 
expected  due  to  the  amount  of  total  energy  available  and  indicates  that  vibrational  energy  is  more 
effective  than  other  types  of  energy  in  promoting  reaction. 

The  works  of  Levine,  et  al.  describe  in  detail  a  compact  method  to  quantify  the 
®5~ancemen^  chemical  reactions  due  to  reactant  internal  energy  using  a  single  parameter,  X. 
When  A,  —  0,  then  there  is  no  vibrational  favoring,  and  the  reaction  probability  is  increased  by  all 
types  of  energy  equally.  This  is  known  as  the  "prior"  case  and  reflects  the  statistically  most  likely 
situation  in  the  absence  of  any  constraint  and,  therefore,  the  maximum  entropy  case  A 
discussion  of  this  approach,  as  applied  to  DSMC,  may  be  found  in  Marriott*  and  Gallis.6  A  large 
negative  value  of  X  indicates  that  the  reaction  is  strongly  vibrationally  favored.  It  has  been  shown 
that  many  ^exchange  reactions,  especially  of  hydrogen  and  halogen-bearing  species,  are  highly 
vibrationally  favored,  while  some  others  show  no  such  favoring. 

For  air  chemistry,  the  two  Zel'dovich  atom  exchange  reactions  are  important: 

N2(v)  +  O  <-»  NO  +  N  (1) 

NO(v)  +  0o02+N  (2) 

For  each  of  these  reactions  (exothermic  direction,  <-),  a  measurement  of  vibrational 
energy  disposal  to  the  products  is  available7*8  at  room  temperature,  in  both  cases  showing  a 
nearly  statistical  distribution.  Detailed  balance  then  dictates  that  in  the  endothermic  direction 
(->)  these  two  reactions  are  not  (or  are  weakly)  vibrationally  favored.  Recent  quasi-classical 
trajectories  (QCT)  calculations9  of  state-to-state  reaction  cross  sections  for  Reaction  (1)  provide 
valuable  detailed  information,  yet  also  confirm  this  qualitative  conclusion.  The  QCT  results 
show  that,  for  a  given  translational  and  rotational  temperature,  the  reduction  in  reaction  rate 
constant  due  to  decreasing  vibrational  temperature  of  the  reactants  is  about  what  would  be 
expected  for  a  case  with  a  slight  vibrational  favoring  of  about  X  =  -2.  Thus,  traditional  rarefied 
flow  data  are  not  the  only  source  for  validation  of  nonequilibrium  chemistry  models.  In  die  case 
ot  the  Zeldovich  reactions,  literature  results  on  the  vibrational  excitation  of  products  for  the 
reverse  reactions  provide  a  valuable  test  of  the  degree  of  vibrational  favoring  of  the  forward 
reactions. 


DISSOCIATION 

Collision-induced  dissociation  reactions  of  diatomic  molecules  are  often  considered  to  be 
vibrationally  favored,  yet  unambiguous  information  on  the  degree  of  favoring  is  scant  Since 
dissociation  reactions  under  conditions  of  significant  vibrational  nonequilibrium  can  be  the 
controlling  chemical  step  in  the  systems  mentioned  here  in  the  introduction,  an  ability  to  model 
these  reactions  with  the  correct  amount  of  vibrational  favoring  is  needed. 
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The  state-specific  cross  sections  for  the  dissociation  of  H2(v)  colliding  with  argon  have 
been  calculated11  using  QCT.  For  this  case  the  vibrational  favoring  for  this  reaction  is  quite 
strong,  corresponding  to  X  of  approximately  -12,  but  H2,  with  its  very  widely  spaced  vibrational 
and  rotational  levels,  is  not  necessarily  a  good  model  for  other  diatomic  species.  A  molecular 
beam  study  of  H2+(v)  dissociating  on  collision  with  He  also  showed  a  strong  vibrational 
favoring12,  reflecting  a  X  of  about  -8.  Several  state-specific  cross  sections  for  the  dissociation  of 
H2  colliding  with  H  have  been  calculated  using  QCT; 1:1  the  results  demonstrate  a  strong 
enhancement  of  dissociation  probability  by  internal  energy,  especially  vibrational  energy. 
Unfortunately,  neither  state-specific  dissociation  measurements  nor  QCT  calculations  are 
available  for  02  and  N2.  Instead,  some  shock  tube  data  on  the  nonequilibrium  dissociation  rate 
coefficient  and  incubation  time  have  been  obtained  for  these  species,  and  a  number  of  studies 
have  attempted  to  deduce  information  on  the  degree  of  vibrational  favoring  based  on  the  these 
macroscopic  quantities. 

It  has  long  been  recognized  that  vibrational  relaxation  and  dissociation  processes  are  fully 
coupled  in  shocks  and  other  rarefied  gas  flows.  Even  without  any  vibrational  favoring, 
dissociation  will  always  occur  more  rapidly  from  molecules  at  higher  vibrational  levels 
compared  with  molecules  at  lower  vibrational  levels.  This  is  because  of  the  much  greater  amount 
of  total  energy  available  in  the  collision  if  the  two  molecules  have  the  same  translation/rotation 
temperature  or  energy.  So,  the  dissociation  process  is  constantly  removing  high  v  (vibrational 
level)  molecules  and  driving  the  gas  away  from  a  Boltzmann  distribution  of  vibrational  energy, 
while  vibration-to-translation  (YT)  relaxation  is  constantly  reshuffling  the  vibrational 
distribution  function  (VDF)  toward  a  Boltzmann  distribution.  The  competition  between  these 
processes  is  known  as  coupled  vibration  dissociation  (CVD)  and  is  what  determines  the  kj*,  the 
quasi-steady  nonequilibrium  dissociation  rate  constant  that  is  observed  in  shocks  and  is  reduced 
from  the  equilibrium  value  of  k^.  An  approximate  value  of  pxd,  where  p  is  pressure  and  is  sn 
approximate  characteristic  time  for  dissociation,  is  included  in  Figure  1,  which  illustrates  that, 
for  02-argon,  CVD  is  important  at  temperatures  above  about  4000  K,  where  t<$  and  xv  (the 
vibrational  relaxation  time)  become  of  comparable  magnitude.  Therefore,  any  model  of  a  shock 
tube  dissociation  measurement  must  include  information  on,  or  assumptions  about,  both  the 
state-specific  dissociation  cross  sections  and  the  state-specific  vibrational  relaxation  cross 
sections.  A  number  of  studies  have  addressed  this  problem  and  many  have  modeled  shock  tube 
dissociation  data  using  various  dissociation  and  vibrational  relaxation  models. 14,18-25 

VIBRATIONAL  RELAXATION  SIMULATIONS 

The  case  chosen  for  these  studies  is  02  in  an  isothermal  heat  bath  of  argon.  This  has  been 
the  subject  of  a  number  of  previous  modeling  studies,  ^-23  primarily  because  of  its  simplicity 
and  the  availability  of  experimental  shock  tube  measurements  of  k/  and  incubation  time  for  this 
case  over  a  range  of  temperatures  by  Wray26  and  Camac  and  Vaughan.27  The  molecules  have  the 
translational  and  rotational  degrees  of  freedom  in  equilibrium  at  some  temperature  T,  but  at  time 
t  =  0,  the  molecules  have  a  vibrational  temperature  (Boltzmann  distribution)  of  200  K.  The  02  is 
considered  to  be  very  dilute  so  that  the  bath  temperature,  T,  remains  constant  and  02-02 
collisions  may  be  neglected.  In  this  way,  studies  of  the  effect  of  the  key  features  of  CVD  can  be 
performed  in  a  very  simplified  context.  All  simulations  presented  here  have  used  T(bath)  =  7000 


The  code  used  here  follows  the  basic  methodology  of  Bird.211  The  relaxation  logic  is 
based  on  the  Particle-Selection  Prohibiting  Multiple  Relaxation  method  of  Haas,  et  al.29  More 
complete  discussion  of  the  current  code  will  be  given  in  Wadsworth’s  “Comparison  of 
Vibration-Translation  Models  for  DSMC  Method.”1  For  the  02-argon  isothermal  bath  of  interest, 
one  needs  only  to  consider  collision  events  involving  one  02  molecule  with  one  argon  atom,  with 

+  Wadsworth,  D.C.,  Raytheon  STX,  AFRL,  Edwards  AFB,  CA,  In  preparation. 
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only  02  dissociation  and  02  vibrational  relaxation  events  being  allowed.  For  consistency 
collision  parameters  for  02  with  argon  are  based  on  the  variable  soft  sphere  (VSS)  values  given 
by  Koura.  The  02  vibration  levels  are  based  on  the  bounded,  quantized,  simple  harmonic 
oscillator  mode!,  with  a  vibrational  temperature  ©v  =  2274  K,  yielding  27  bound  vibrational 
levels.  The  various  models  of  dissociation  and  VT  exchange  are  implemented  into  the  basic  code 
with  only  mmor  changes  to  the  collision  logic. 


0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.1 

j- 1/3 

Figure  1 

Predicted  Correlation  of  pxv  with  Temperature  for  02-argon  from  Various  Sources 

An  approximate  value  ofpxd  (average  time  to  dissociation)  is  included  for  comparison 


It  is  worth  stating  that  the  only  calibration  used  for  many  VT  models  is  that  of  matching 
the  overall  vibrational  relaxation  time  (xv)  of  the  gas  at  a  given  temperature  or  range  of 
temperatures  to  that  predicted  by  a  measured  correlation.  Three  such  correlations  for  02-argon 
collisions  found  in  the  literature  are  shown  in  Figure  1.  The  correlation  of  Camac16  is  based  on 
data  which  extend  from  1200  to  5000  K.  The  correlation  of  Millikan-White1*  (MW)  is  based  on 
data  from  White  and  MillikanS 1  from  1000  to  3000  K  as  well  as  the  data  of  Camac  The 
coirelafron  of  Losev1 7  is  based  on  data  from  4000  to  10,000  K  as  well  as  previous  data 
Unfortunately  the  validity  of  extrapolating  any  correlation  to  higher  temperatures  is 
questionable.  In  addition,  the  assumptions  of  the  Landau-Teller  model  (harmonic  oscillator,  only 
Av  -  1,  Boltzmann  vibrational  distribution),  which  are  used  to  reduce  shock  data  to  obtain  a 

V  *°r  Tv’  k^_gm  t0  krefk  down  at  T  >  2000  K.  Thus,  a  correct  value  for  xv  is  problematic  and 
will  be  seen  to  have  an  effect  on  kd*. 

r»eiv>rr<  ^  vario1us  m®thJ°,ds  °f  modeling  VT  transfer  that  have  been  incorporated  into  the 
DSMC  code  are  described  here  bnefly.  Model  I  is  the  Borgnakke-Larsen  (BL)  method, 32  which 
simulates  VT  relaxation  by  allowing  a  change  in  vibrational  state  (or  energy)  every  Z  collisions. 
Here,  Z  is  the  vibrational  collision  number,  based  on  the  translation-rotation  temperature  of  the 
g  ,,.s?  overaft  relaxation  time  Xy  matches  the  "known"  experimental  value.  When  a 

collision  is  permitted  to  undergo  VT  relaxation,  the  post-collision  state  (or  energy)  of  the  02 


5 


molecule  is  sampled  from  an  equilibrium  distribution  based  on  equipartition  of  the  sum  of 
translational  and  vibrational  energy  of  the  collision  among  the  available  modes. 

In  contrast  to  the  BL  model,  one  may  dispense  with  Z  and  simulate  VT  transfer  using 
state-to-state  transition  probabilities.  This  was  first  implemented  for  DSMC  by  Boyd33  and  later 
by  Abe,34  who  formulated  probabilities  for  Av  =  1  transitions  as  a  function  of  collision  velocity 
and  initial  vibrational  level  in  such  as  way  as  to  reproduce  the  correlation  between  xy  and  gas 
temperature.  It  seems  that,  as  the  gas  temperature  increases,  the  velocity-dependent  transition 
probabilities  in  this  model  do  not  produce  equipartition  of  energy,  due  to  approximations  made  in 
the  derivations  (see  Choquet35).  However,  for  the  case  of  a  thermal  bath,  the  velocity-dependent 
probability  is  not  really  needed.  Abe34  has  presented  a  model  with  constant  transition  probability 
for  each  vibrational  level  (independent  of  collision  energy).  (Note  Abe’s  Equation  32  in 
Reference  34  is  valid  for  hard  sphere  only).  This  will  correctly  equipartition  energy  and  can  be 
set  to  produce  the  desired  value  of  xv,  except  at  very  high  temperatures,  where  probabilities  for 
high  v  states  exceed  unity  and  where  multi-quantum  transitions  are  important.  This  model  has 
been  labeled  as  Model  II. 


A  third  model  for  VT  transfer  used  in  this  study  is  that  first  presented  by  Koura^L3^  for 
Monte  Carlo  implementation.  The  probability  equations  are  given  in  Heidrich,  et  al.,37  who 
referred  to  it  as  Improvement  to  Forced  Oscillator,  Impulsive  Transfer  Semiclassical  (ITFITS).  It 
is  based  on  the  forced  harmonic  oscillator  approximation  (FHO),  as  originally  derived  by 
Kemer38  and  Treanor,3^  and  includes  multi-quantum  VT  transfer.  Koura  implements  the  cross 
section  (as  a  function  of  relative  velocity  and  initial  vibrational  level)  for  VT  transfer  for  a 
diatom  colliding  with  an  atom  in  his  Monte  Carlo  simulations  of  02-argon.  The  same  expression 
(in  a  slightly  different  form)  for  VT  probabilities  was  presented  by  Adamovich,  et  a/.40  who  use 
the  probabilities  to  compute  state-specific  rate  coefficients  for  vibrational  transfer.  The  work  by 
Adamovich,  et  al.  also  presents  the  formulas  for  probability  of  vibration- vibration  transfer  and 
the  resulting  rate  coefficients.  It  includes  an  extensive  validation  of  the  FHO  model  by 
comparing  the  results  for  N,-N2  and  02-02  relaxation  rate  constants  with  semi-classical  3-D 
trajectory  results  by  Billing4  H42  over  a  range  of  temperatures  and  vibrational  levels.  In  addition, 
the  value  of  Xy  predicted  by  FHO  for  N2-  N2  and  02-argon  is  shown  to  be  in  reasonable 
agreement  with  measured  values  over  a  range  of  temperatures. 

The  FHO  or  ITFITS  model  has  two  adjustable  parameters,  the  steepness  parameter  of  the 
potential,  a,  and  a  steric  factor  (STF).  For  02-argon,  Koura  estimates  the  values  of  a  and  STF 
that  result  in  exact  agreement  with  the  MW  estimation  of  xv  in  the  range  300-2000  K,  but 
diverges  at  T  >  2000  K.  The  ITFITS  model  for  VT  transfer  has  been  implemented  in  the  DSMC 
code  using  the  value  recommended  by  Koura  for  a  and  STF  and  will  be  called  Model  III. 

Some  results  of  heat  bath  calculations  for  02-argon  using  the  three  VT  models  are 
presented  in  Figures  2  and  3.  The  mean  vibrational  energy  of  a  gas  may  be  defined  as: 


£evn. 


V 


(3) 


where  Ev  is  the  energy  of  each  vibrational  level  (assuming  that  the  ground  level  has  Ev  =  0),  and 
Nv  is  the  number  of  molecules  that  populate  that  vibrational  level.  Figure  2  shows  several 
examples  of  (v)  versus  time  in  the  heat  bath,  demonstrating  that,  in  the  absence  of  dissociation, 
all  VT  models  used  in  this  study  achieve  the  correct  equilibrium  value  of  (v),  but  vary  in  the 
length  of  the  relaxation  time.  The  shape  of  all  the  curves  in  Figure  2  follows  that  predicted  (not 
shown)  by  the  Landau-Teller  (LT)  equation:  d(v)/dt  =  ((v)eq-(v»/xv  (note  that  this  equation 
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provides  the  definition  of  xv).  Here,  <v)eq  is  the  equilibrium  mean  vibrational  energy  at  the  heat 
bath  temperature.  Since  die  various  VT  models  reproduce  the  shape  of  the  LT  equation  the 
overall  vibrational  relaxation  predicted  by  any  model  may  be  approximated  by  some  value  of  Tv 
that  causes  the  LT  equation  to  match  the  plot  of  (v)  versus  time.  Figure  2  demonstrates  the  range 
of  uncertainty  m  xv  at  7000  K.  6 


Figure  2 

Relaxation  of  the  Mean  Vibrational  Energy  in  the  Isothermal  Bath 

The  results  using  Model  I  are  not  shown,  since  for  each  value  of  xv,  the  curve  for  Model  I  lies  on 

top  of  the  curve  for  Model  II 


When  Models  I  and  II  use  the  same  value  of  tv,  the  relaxation  of  <v>  versus  time  is  the 
same.  Yet,  at  short  times,  the  actual  VDFs  produced  by  Models  I  and  II  can  differ  substantially 
Figure  3  demonstrates  a  simulation  using  Model  I  and  Model  II  with  the  same  value  of  x  Figure 
3(a)  shows  tiie  resulting  VDFs  at  an  early  time  during  the  vibrational  relaxation,  t »  0.1  *x  .  (The 
MW  value  for  xv  is  used  m  this  case).  Assuming  that  Model  II  is  closer  to  reality,  the  result  of 
Model  I  shows  a  systematic  bias  in  VDF  in  the  short  time  (few  collision)  regime.  One  can  see  by 
wrth hie  showing  a  Boltzmann  VDF  for  7000  K how  Model  I  creates  a  bimodi 

Y  . \ rrvr0^^00} ,10ns  dIt  allowed  to  undergo  vibrational  relaxation  immediately  jump  to 
the  VDF  of  the  bath.  J  * 
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One  may  define  a  “vibrational  temperature,”  Tv,  by: 

£Evexp(-Ev/kTv) 

<v>  =  -^= - — 

£exp(-Ev/kTv) 

V 

(where  k  is  Boltzmann’s  constant),  and  solving  for  Tv  based  on  the  value  of  at  <v)  any  point  in 
time.  Figure  3(b)  shows  that,  based  on  this  definition,  the  T  produced  by  the  two  models  follows 
the  same  progression  in  time.  However,  if  Tv  is  defined  as: 


T  = - -  (5 

V  ln(N0/N,)  P- 

as  is  often  the  case  (where  N0  is  the  population  in  v  =  0  and  N,  is  the  population  in  v  =  1),  Figure 
3(c)  shows  that  Model  I  produces  a  different  value  of  Tv  compared  with  Model  II  at  short  times, 
due  to  the  differences  in  VDF.  It  must  be  emphasized  at  this  point  that  a  gas  whose  VDF  does 
not  follow  a  Boltzmann  distribution  does  not  have  a  well-defined  vibrational  temperature,  and 
any  presentation  of  an  approximate  Tv  of  such  a  gas  is  not  meaningful  without  an  accompanying 
statement  of  how  the  approximate  Tv  is  defined.  The  systematic  bias  of  the  BL  model 
complicates  any  attempt  at  quantitative  comparisons  between  results  in  the  literature  of  Tv. 


vibrational  level  (v) 

Figure  3(a) 

Vibrational  Distribution  Function  Predicted  by  Model  I  and 
Model  II  at  t  =  0.1  *xv  (using  tv  MW) 
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e<*n-  <5>  Tv  .  eqn.  (4) 


figure  3(b) 

Vibrational  Temperature  as  a  Function  of  Time  using  Equation  (4) 
for  Model  I  and  Model  II 


.  Figure  3(c) 

Vibrational  Temperature  as  a  Function  of  Time  using  Equation 
for  Model  I  and  Model  II  V 
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INTRODUCTION  OF  DSMC  DISSOCIATION  MODELS 

As  discussed  in  the  previous  section,  many  methods  for  modeling  dissociation  of 
diatomics  have  been  presented  in  the  past.  The  number  of  methods  that  have  been  used  in  the 
context  of  DSMC  simulations  is  more  limited.  The  method  introduced  by  Bird43  to  simulate 
chemical  reactions  in  DSMC  calculations  is  widely  utilized.  It  expresses  the  probability  of 
dissociation  as  a  function  of  the  total  collision  energy  (TCE)  and  uses  adjustable  parameters  to 
reproduce  measured  dissociation  rate  coefficients  as  a  function  of  temperature  under  equilibrium 
conditions.  Recently,  Bird  introduced  a  different  DSMC  method  for  chemical  reactions,28’44’45 
known  as  the  exact  available  energy  (EAE)  approach.  It  performs  a  BL  vibrational  redistribution 
on  every  collision  that  has  sufficient  energy  to  dissociate  and,  using  “notional”  vibrational  levels 
that  extend  above  the  dissociation  energy  (D),  chooses  a  post-collision  level  based  on  the  total 
pre-collision  energy.  If  the  post-collision  level  is  above  D,  then  a  dissociation  reaction  is 
assumed  to  have  taken  place.  Under  equilibrium  conditions,  this  method  has  been  shown  to 
reproduce  fairly  closely  the  measured  dissociation  rates  for  02-02  and  N2-N2  as  a  function  of 
temperature.  Like  the  TCE  method,  the  EAE  is  based  on  the  total  energy  of  a  collision,  with  no 
favoring  of  vibrational  energy.  The  EAE  dissociation  model  has  been  implemented  in  die  DSMC 
code.  More  details  are  provided  in  Reference  46.  Figure  4  shows  the  equilibrium  dissociation  rate 
coefficient  predicted  by  the  EAE  model  for  02-argon  collisions. 

Haas  and  Boyd20  introduced  the  vibrationally  favored  dissociation  (VFD)  model  which 
modifies  the  TCE  model  and  allows  for  favoring  of  vibrational  energy  in  the  dissociation 
probability  through  an  adjustable  (empirical)  parameter,  <j».  Hash  and  Hassan47  have  proposed  a 
modified  version  of  the  VFD  model  where  the  value  of  <|>  is  determined  without  experimental 
input.  Ivanov,  et  a/.,48  have  implemented  a  dissociation  model  in  their  DSMC  code  that  uses 
probabilities  based  on  state-specific  rate  constants  (as  formulated  by  Wamatz49)  that  include  a 
certain  (fixed)  amount  of  vibrational  favoring.  Macheret  and  Rich18  have  presented  a  threshold 
line  model  to  obtain  dissociation  rate  coefficients  for  given  (unequal)  translation/rotation  and 
vibrational  temperatures.  Recently,  Boyd50  has  adapted  this  model  for  a  DSMC  implementation, 
using  one  probability  expression  for  dissociation  from  high  vibrational  levels  and  a  different 
(reduced)  probability  expression  for  low  vibrational  levels.  The  reduction  in  probability  for  lower 
levels  means  that  this  model  is  vibrationally  favored.  The  probabilities  include  a  single 
adjustable  STF  for  calibrating  the  results  for  equilibrium  against  measured  rate  constants. 

An  approach  to  modeling  vibrationally  favored  dissociation  based  on  the  formalism  of 
Levine4  (as  discussed  above  for  exchange  reactions)  was  used  by  Kiefer,19  Gonzales,22*23  and 
Shizgal24  for  continuum  simulations  of  nonequilibrium  gases.  For  particle  simulations, 
Marriott,5  Gallis,0  and  Koura21  have  presented  adaptations  of  the  Levine  approach;  all  of  these 
include  an  adjustable  parameter,  X,  that  indicates  the  degree  of  vibrational  favoring,  but  is 
formulated  in  different  ways.  The  concept  for  a  Monte  Carlo  implementation  is  to  have  a  prior 
state-specific  cross  section,  a0,  that  reflects  the  probability  for  dissociation  from  a  given  level  as 
a  function  of  translational  energy.  The  prior  cross  section  is  intended  to  express  the  probability 
with  no  constraint  (such  as  vibrational  favoring)  and  thus  the  condition  of  "maximum  entropy." 
Then  the  prior  cross  section  is  multiplied  by  a  factor  exp(-A.f(Ev,  Etot))  to  simulate  the  constraint 
of  vibrational  favoring,  where  f(v)  is  related  to  the  fraction  of  collision  energy  that  is  in  the 
vibrational  mode.  (Note  that  there  is  a  sign  change  between  Levine's  papers  and  some  of  the 
more  recent  formulations.  Originally,  the  formulas  were  such  that  a  negative  value  of  X 
represented  vibrational  favoring.  Some  recent  versions,  including  Koura's  and  References  19  and 
22-24  have  changed  the  sign  such  that  a  positive  value  of  X  represents  vibrational  favoring.) 
Unfortunately,  the  prior  cross  section  for  a  dissociation  reaction  (which  in  principle  comes  from 
counting  possible  post-reaction  states)  is  dependent  on  the  intermolecular  potential  that  is 
assumed53  and  a  general  expression  is  not  known.  Kafri  and  Levine54  have  suggested  one 
approximation  as  ct°  =  A(E-D)2/(Et)1/2. 
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Figure  4 

Rate  Constants  for  Dissociation  of  02  by  Argon 
The  DSMC  results  represent  true  equilibrium  rate  constants  (the  VDF  of  the  molecules  in  the 
simulation  is  constrained  to  Tv  =  T).  The  Marrone  value  (Equation  7),  as  explained  in  the  text,  is 
estimated  for  true  equilibrium,  and  the  Gupta  value  (Equation  8)  and  Baulch  value  (Equation  9) 
are  based  on  experimental  measurements  that  reflect  some  nonequilibrium 


Although  the  continuum  versions  mentioned  above  (References  19,  22-24)  make  use  of 
the  Kafri  and  Levine54  suggestion  for  the  prior  cross  section,  which  includes  the  total  collision 
energy,  Koura^1  uses  a  different  expression,  ct°  =  ael(l-(D-Ev)/s).  Here  <yel  is  the  elastic 
collision  cross  section  based  on  the  symmetrized  relative  collision  velocity,  Ey  is  the  vibrational 
energy,  and  e  is  the  relative  translational  energy.  This  prior  cross  section  does  not  include  any 
contribution  to  dissociation  probability  from  rotational  energy,  which  may  not  be  physically 
realistic  for  some  cases.  The  full  expression  for  dissociation  cross  section  used  by  Koura  is: 

crd  =  cjel(l-(D-Ev)/s)STF  exp[X(Ev/D  - 1)] 

for  s>(D-Ev)  (6) 

and  ad  =  0  otherwise.  Koura  recommends  STF  =  1  and  X  =  2  for  02-argon.  In  the  DSMC  code 
used  for  the  present  parametric  studies,  this  expression  has  been  implemented.  The  STF  was  set 
to  011j/or  Pri°r  case  =  0)-  When  the  value  of  X  is  changed,  however,  the  value  of  the 
adjusted  accordingly  so  that  the  degree  of  vibrational  favoring  does  not  change  the 
equilibrium  dissociation  rate  constant.  For  instance,  when  X  =  2,  the  STF  =  2.4.  Since  Koura 
always  uses  STF  =  1 ,  his  k<je£l  will  be  lower  than  this  one  for  X  —  2.  As  shown  in  Figure  4,  under 
equilibrium  calculations  the  rate  constant  predicted  by  the  Koura  model  for  02~urgon  and  that 
from  the  EAE  model  are  virtually  identical,  but,  as  seen  below,  the  rate  constant  predicted  under 
nonequilibnum  (shock  tube)  conditions  by  the  Koura  model,  even  for  the  X  =  0  (prior)  case,  is 
lower  than  that  of  the  EAE  model.  This  may  be  because  Koura's  expression  for  the  prior  cross 
section  is  not  truly  representative  of  the  maximum  entropy  case. 
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DISSOCIATION  SIMULATION  RESULTS 

Ideally,  one  would  like  to  validate  a  nonequilibrium  dissociation  model  in  the  following 
steps.  First,  ensure  that  the  model  predicts  the  correct  rate  constant  under  fully  equilibrium 
conditions.  Then,  compare  the  model  predictions  for  the  quasi-steady  (nonequilibrium) 
dissociation  rate  constant  and  incubation  time  in  a  shock  with  measured  values.  Unfortunately 
both  of  these  steps  are  far  from  straightforward.  The  various  literature  measurements  of  the 
dissociation  rate  constant  for  a  given  system  (02-argon  is  used  here)  show  a  considerable  amount 
of  scatter,  which  tends  to  increase  at  higher  temperatures.  In  addition,  most  of  these  measured 
values  are  from  gases  with  varying  amounts  of  nonequilibrium,  and  thus  do  not  truly  measure  the 
equilibrium  rate  constant.  Marrone  and  Treanor^  have  stated  that  vibrational  nonequilibrium 
introduces  the  pre-exponential  temperature  dependence  in  the  Arrhenius  expression  for  the 
dissociation  rate  constant,  and  recommended 

kdeq  =  9e  1 4  *  exp(-5 93 8 0/T)  cm3/mole*s  (7) 

as  the  equilibrium  rate  constant  expression  for  02-argon.  This  expression  agrees  with  the  shock 
tube  observations  of  Camac  and  Vaughan2?  up  to  4000  K,  but  at  higher  temperatures  diverges  as 
the  shock  tube  data  reflect  the  effect  of  nonequilibrium  as  the  characteristic  time  for  dissociation 
becomes  comparable  to  xv.  Kondratiev  and  Nikitin55  have  stated  that  the  rate  coefficient  for  CVD 
will  be  proportional  to  exp(-D/kT)  in  the  low  temperature  limit  and  proportional  to  T'exp(- 
D/kT)  m  the  high  temperature  limit.  Haas  and  Boyd20  have  adopted  the  above  expression 
recommended  by  Marrone  to  calibrate  their  VFD  model  for  the  equilibrium  case. 

It  should  be  noted  that  the  two  dissociation  models  employed  in  this  study  do  not  use 
adjustable  parameters  to  set  the  equilibrium  rate  constant,  although  an  STF  could  be  added  to 
reduce  the  predicted  rate  constant,  if  desired.  It  happens  that  the  two  models  produce  nearly 
identical  equilibrium  rate  constants  for  02-argon,  as  shown  in  Figure  4.  Both  are  a  factor  of  2-3 
higher  than  the  Mamme  recommendation,  and  up  to  a  factor  of  10  higher  than  the  expression 
recommended  by  Gupta5 1 : 

kj  =  3.61el8*T'°exp(-59400/T)  cm3/mole*s  (8) 


or  Baulch52: 

k„  =  1 . 8  Oe  1 8  *  T!  °exp(- 5938 0/T)  cm3/mole*s  (9) 

Both  the  Gupta  and  Baulch  recommendations  are  based  on  all  available  measurements  and 
presumably  reflect  the  effects  of  nonequilibrium  at  temperatures  above  4000  K.  However  the 
Baulch  recommendation  is  a  factor  of  2  lower  than  Marrone's14  in  the  3000-4000  K  range  Also 
more  recent  shock  tube  data  by  Jerig  et  aL*  fit  the  C>2-argon  dissociation  rate  constant  for  the 
range  2400  to  4100  K  to: 

kj  =  1.6el8*T10exp(-59380/T)  cm3/mole*s  (10) 

identical  to  the  Baulch  recommendation.  Thus,  the  proper  value  for  kd'q  to  match  the 
DSMC  model  to  is  ambiguous. 

The  second  step  of  the  validation  -  comparing  nonequilibrium  rate  constants  with  shock 
measurements,  is  complicated  by  the  influence  of  vibrational  relaxation,  as  previously  discussed, 
rigures  5(a)  and  5(b)  show  parametric  studies  to  illustrate  the  effects  of  various  assumptions  on 
the  results  of  shock  tube  simulations.  The  number  of  02  molecules  in  the  simulation  cell  is 
shown  as  a  function  of  time  in  a  semilog  plot.  The  slope  of  the  linear  (quasi-steady)  portion 
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yields  the  value  of  kd  *  for  each  simulation.  The  intercept  of  the  extrapolated  line  with  a 
honzontal  line  at  n(02)  =  1  gives  the  value  of  Tinc,  (the  incubation  time).  Figure  5(a)  shows  how 
V  model  affects  die  simulation  result  for  kd  *.  Even  if  the  dissociation  model  is  identical 
(the  EAE  model  is  used  for  both  cases  shown  in  the  figure),  using  a  VT  model  (Model  II  is  used 
here)  that  has  a  different  effective  tv  will  result  in  a  difference  inkd  *. 


Figure  5(a) 

Effect  of  Ty  on  Rate  of  Nonequilibrium  Dissociation 


Figure  5(b)  shows  the  effect  of  vibrational  favoring  in  the  dissociation  model.  The  four 
cases  shown  all  use  Model  III  for  VT,  but  use  EAE28  dissociation  or  Koura28  dissociation 
with  A  0,  1,  and  2  to  vary  the  amount  of  vibrational  favoring,  which  results  in  a  change  in  k,* 
°*  —PH*  a  ®ctor.°f  Also  included  in  the  figure  for  comparison  is  an  indication  of  the 
equilibrium  dissociation  rate  constant  predicted  by  the  current  models  to  re-emphasize  that  Jr*  is 
significantly  reduced  from  kdeq  even  with  no  vibrational  favoring  (EAE  model).  Table  1 
summarizes  the  results  for  kd*  and  Tinc  for  a  number  of  simulations.  Also  included  are 
expenmentaUy  measured  values  for  Tinc  and  k/,  but,  since  dissociation  models  typically  have  an 
adjustable  stenc  factor,  only  the  ratio  of  kdeq/kd*  is  truly  significant,  and  a  measured  value  of  k,* 
without  a  known  value  of  kd  cannot  be  used  to  distinguish  among  various  models. 


time  (s) 

Figure  5(b) 

Effect  of  Varying  Amount  of  Vibrational  Favoring  on  Nonequilibrium  Rate  of  Dissociation 

A  line  showing  the  rate  of  equilibrium  dissociation  is  included  for  comparison 


Table  1.  Simulation  Results  for  Nonequilibrium  Dissociation  of  P2  in  Argon  at  7000  K 


*inc  (^sec) 

ka*  (xlO19  m3/s) 

kdeq  /kd* 

Model  II,  MW  xv,  EAE  dissociation28 

1.4 

5.2 

1.5 

Model  II,  Losev  Xy,  EAE  dissociation28 

2.2 

4.3 

1.9 

Model  III,  EAE  dissociation28 

2.1 

4.4 

1.8 

Model  III,  Koura21  dissociation,  X  =  0 

3.0 

2.8 

2.8 

Model  III,  Koura21  dissociation,  X  =  1 

3.8 

2.0 

3.9 

Model  III,  Koura21  dissociation,  X  =  2 

4.5 

1.6 

5.0 

Measured  -  Wray26 

4-5 

8 

N/A 

Figure  6  shows  the  quasi-steady  state  VDF  for  the  same  simulations  presented  in  Figure 
5(b).  For  each  vibrational  level,  f(v)/fq(v)  is  plotted  -  the  fraction  of  02  molecules  in  that  level 
divided  by  the  fraction  that  would  be  present  in  a  Boltzmann  distribution  at  7000  K.  Notice  that, 
because  fi (vyf^v)  dips  below  1  for  the  middle  and  high  vibrational  levels,  it  must  be  greater  than 
1  for  the  lowest  levels.  This  behavior  is  predicted  by  Haug,  et  al.  for  the  H2-argon  system.57  This 
figure  may  be  compared  with  Figure  7  of  Reference  20  which  plots  the  same  quantity  for  the 
TCE  and  VFD  dissociation  models  for  02-argon,  also  at  7000  K  (Reference  20  uses  anharmonic 
oscillator  energy  levels,  but  the  behavior  at  the  same  value  of  Ey/D  should  be  the  same).  The 
TCE  result  from  that  plot  and  the  EAE  result  in  Figure  6  here  should  be  the  same,  except  for  the 
following  two  differences.  First,  the  overall  xy  for  the  simulations  shown  in  Figure  6  (ITFITS  VT 
model)  is  significantly  greater  than  the  xy  value  (based  on  Camac’s  correlation)  assumed  in 
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Reference  20.  Second,  the  dissociation  models  in  Reference  15  have  been  adjusted  to  agree  with 
Marrone's  recommendation  for  kdeq,  which  is  a  factor  of  2.5  lower  than  produced  by  the  EAE 
model.  Both  of  these  differences  will  tend  to  drive  the  quasi-steady  state  VDFs  shown  in  Figure 
6  farther  away  from  a  Boltzmann  distribution  than  the  cases  shown  in  Reference  20.  This 
provides  another  example  of  the  sensitivity  of  nonequilibrium  simulation  results  to  die  details  of 
the  model  inputs.  The  trend  that  can  be  seen  in  Figure  6  is  that  increasing  the  amount  of 
vibrational  favoring  in  the  dissociation  model  will  drive  the  quasi-steady  state  VDF  down  for  the 
very  highest  vibrational  levels,  but  up  in  the  intermediate  levels. 


vibrational  level 


Figure  6 

Quasi-steady-state  Vibrational  Distribution  Function  with  Different  Amounts  of 

Vibrational  Favoring 


CONCLUSIONS 

For  the  atom  exchange  reactions  relevant  to  air  chemistry  (Zel'dovich  reactions),  there 
exists  reasonable  validation  that  a  chemistry  model  based  on  the  total  collision  energy  without 
vibrational  favoring  is  adequate.  For  dissociation  reactions,  a  number  of  studies  have  used  shock 
tube  data  of  nonequilibrium  dissociation  rate  constants  to  show  that  a  vibrational  favoring  effect 
is  very  likely.  More  detailed  model  validation  cannot  be  obtained  from  shock  data,  however 
without  additional  confidence  in  the  estimated  value  of  k^  and  in  the  VT  model,  since  the 
uncertainties  in  these  quantities  cause  effects  of  about  the  same  magnitude  as  the’  effects  of 
variations  in  the  dissociation  model.  It  seems  likely  that  a  fairly  simple  computational  model  will 
be  able  to  capture  the  important  aspects  of  chemically  reacting  nonequilibrium  flows.  However 
true  validation  of  the  chemistry  models  used  in  DSMC  will  require  detailed  state-specific 
measurements  or  QCT  or  other  type  of  calculation  of  rate  constants  or,  better  yet,  energy- 
dependent  cross  sections.  Currently  available  data  do  not  permit  an  evaluation  of  which  model 
best  simulates  nonequilibrium  air  dissociation  reactions.  This  study  emphasizes  the  need  for  a 
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careful  evaluation  of  any  proposed  verification  experiment  as  to  whether  it  will  provide  sufficient 
information  to  distinguish  between  different  DSMC  chemistry  models. 
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PART  II:  EXAMINATION  OF  DSMC  CHEMISTRY  MODELS:  ROLE  OF 
VIBRATIONAL  FAVORING 


INTRODUCTION 

An  accurate  model  of  chemical  reactions  is  a  key  factor  in  obtaining  useful  predictions  of 
radiation  and  surface  heating  around  re-entry  vehicles  or  spacecraft  and  in  understanding  plasma 
processing  environments,  etc.  These  types  of  flows,  however,  typically  display  significant 
thermal  and  chemical  nonequilibrium.  The  degree  of  population  of  high  vibrational  levels  is 
often  the  controlling  factor  in  the  probability  of  dissociation  and  endoergic  exchange  reactions. 

Direct  Simulation  Monte  Carlo  (DSMC)  calculations,  which  simulate  individual  gas 
particle  collisions,  are  well  suited  to  the  inclusion  of  physically  accurate  chemistry  models  based 
on  the  translational  and  internal  energy  states  of  each  molecule.  A  minimum  requirement  for  a 
DSMC  chemistry  model  would  be  to  reasonably  reproduce  known  equilibrium  rate  constants.  A 
second  requirement  is  to  reasonably  reflect  the  relative  probabilities  for  reaction  for  molecules  in 
different  internal  energy  states,  thus  allowing  realistic  modeling  of  gases  where  the  distribution 
of  internal  energies  is  far  from  equilibrium.  A  number  of  proposed  chemistry  models  have  been 
demonstrated  to  meet  the  first  requirement.  The  present  study  concentrates  on  questions  relating 
to  the  second  requirement. 

Some  chemical  reactions  are  known  to  have  their  reaction  probability  enhanced  by 
vibrational  excitation  of  the  reactants.  This  characteristic  will  be  referred  to  as  "vibrational 
favoring."  Any  endoergic  reaction  will  be  enhanced  by  vibrational  excitation,  merely  because  the 
amount  of  total  energy  available  in  the  collision  is  increased  (assuming  that  all  types  of  energy 
contribute  to  the  activation  energy).  "Vibrational  favoring"  refers  to  an  increase  in  reaction 
probability  beyond  what  would  be  expected  due  to  the  amount  of  total  energy  available  and 
indicates  that  vibrational  energy  is  more  effective  than  other  types  of  energy  in  promoting 
reaction. 


EXACT  AVAILABLE  ENERGY  (EAE)  MODEL 


The  method  introduced  by  Bird1  to  simulate  chemical  reactions  in  DSMC  calculations, 
sometimes  referred  to  as  total  collision  energy  (TCE)  model,  is  widely  utilized.  It  uses  adjustable 
parameters  to  reproduce  measured  reaction  rate  coefficients  as  a  function  of  temperature  under 
equilibrium  conditions.  The  number  of  energy  modes  of  the  reacting  molecules  contributing  to 
the  reaction  probability  is  an  input  parameter  in  the  TCE  model.  It  can  range  from  the  sum  of  all 
internal  energy  and  translational  energy  to  only  including  translational  energy  modes. 

Bird  and  Carlson  have  introduced  a  different  DSMC  method  for  chemical  reactions,2’3 
known  as  the  exact  available  energy  (EAE)  approach  (it  has  also  been  referred  to  as  a 
vibrationally  linked  chemistry  approach).  It  does  away  with  reaction  rate  constants  as  such  and 
allows  dissociation  or  atom  exchange  whenever  the  post-collision  vibrational  energy  based  on 
Borgnakke-Larsen  (BL)  redistribution  reaches  or  exceeds  die  activation  energy.  Under 
equilibrium  conditions,  this  method  has  been  shown2’3  to  reproduce  the  correct  order  of 
magnitude  for  the  dissociation  rates  of  O2-O2  and  N2-N2  as  well  as  the  Zel'dovich  air  exchange ' 
reaction  rates  as  a  function  of  temperature.  Like  the  TCE  method,  the  EAE  is  based  on  the  total 
energy  of  a  collision,  with  no  favoring  of  vibrational  energy. 

Other  models  generally  use  input  parameters  that  are  adjusted  to  reproduce  "known" 
equilibrium  rate  constants  that  are  in  fact  not  well  known  at  high  temperatures.  Thus,  one 
interesting  feature  of  the  EAE  model  is  that  it  does  not  use  adjustable  input  parameters. 
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However,  there  is  some  variability  in  results  for  dissociation  due  to  the  use  of  "notional" 
vibrational  levels.  These  are  assumed  vibrational  levels  that  extend  above  the  dissociation  energy 
(D)  with  a  suggested  spacing  either  equal  to  2  or  one-half  of, 3  the  spacing  of  the  two  highest 
bound  vibrational  levels.  The  BL  vibrational  redistribution  method  chooses  a  post-collision  level 
based  on  the  total  pre-collision  energy.  If  the  post-collision  level  is  above  D,  then  a  dissociation 
reaction  is  assumed  to  have  taken  place.  In  this  study,  it  has  been  found  that  the  choice  of  the 
energy  spacing  of  the  notional  levels  has  a  direct  effect  on  the  equilibrium  dissociation  rate 
constant  that  is  predicted  by  the  EAE  model,  with  the  rate  constant  decreasing  as  the  spacing  of 
the  notional  levels  increases.  r  & 

r  be  noted  that  the  probabilities  for  recombination  and  reverse  exchange  reactions 

tor  the  EAE  chemistry  model  are  given  in  Reference  2  based  on  equilibrium  collision  theory  and 
calculated  partition  functions.  But,  since  the  forward  dissociation  probability  depends  on  the 
notional  level  spacing,  the  expression  for  the  recombination  probability  must  be  adjusted 
depending  on  the  spacing  used  for  the  notional  levels  in  order  to  correctly  reproduce  the 
equilibrium  fraction  of  dissociation.  Also,  as  pointed  out  in  Reference  2,  the  EAE  model  for 
dissociation  is  nearly  independent  of  collision  partner  (except  for  small  changes  due  to  reduced 
I?1?  6  hjFdu?phere  (VHS)  collision  parameter  co),  and  so  does  not  allow  for  the  so- 
called  third  body  probability  term.  Smce  the  measured  dependence  on  collision  partners  typically 

the  modef  expenmental  uncertainty>  however,  it  is  difficult  to  assess  the  impact  of  this  future  of 

ATOM  EXCHANGE  REACTIONS 

For  air  chemistry,  the  two  Zel'dovich  exchange  reactions  are  important,  one  of  which  is: 

N2(v)  +  0<->  NO  +  N  (1) 

It  is  known4-6  tfiat  some  endoergic  atom  exchange  reactions  display  a  large  vibrational  favoring 
effect,  while  others  display  virtually  no  preference  for  vibrational  energy  as  compared  to  other 
typ^  of  energy  ft  has  been  argued  previously^  that  infoimation  from  the  vibrational  energy 
distribution  in  the  products  of  the  reverse  (exoergic)  reaction,  along  with  detailed  balance 

S?  it  be  ^ed  t0  ^ermine  whether  an  endoergic  exchange  reaction  is  strongly 
vibrationally  favored  or  not.  Based  on  this  approach,  it  was  predicted  that  the  Zel'dovich 
reactions  will  not  have  significant  vibrational  favoring.  For  a  reaction  where  vibrational  favoring 

model)  wffl  beadequate™8^  m°^  *S  based  on  die  total  collision  energy  (such as  the  EAE 

Recent  quasi-classical  trajectory  (QCT)  calculations*.!*)  0f  state-to-state  reaction  cross 
S  .  P5TdC  valuabledetailed  ^formation  for  Reaction  (1).  In  light  of  this  new  information, 
model  tiff.  tWP5P°fw  a  generalized  collision  energy  (GCE)  model  that  modifies  the  TCE 
adjustable  terms  to  separately  weight  the  amount  of  translational,  rotational, 
and  vibrational  energy  of  the  reactants.  One  measure  of  the  performance  of  the  GCE  model  is 

ArDSMr°o3fS?Ae  diStTbw0n  °f  ^  fo^ogen  molecules  that  are  selected  for  reaction  by 

predicted  by  the  QCT  results  for  two  different  temperature 

for  pr°portlonal  tbe  product  ofthe  relative  reaction  rate  constant 

for  each  vibrational  level  and  the  relative  population  of  that  level).  Excellent  agreement  is  seen. 

a^ionAthat  a  chemistry  model  based  on  the  total  collision  energy 
(  uch  as  the  EAE  model)  will  be  adequate  for  this  reaction,  the  same  two  temperature  conditions 
for  Reaction  (1)  were  simulated  using  a  DSMC  code  with  the  EAE  chemistrymodel.  The  results 
are  shown  in  Figure  1  for  the  equilibrium  case  with  T  =  5000  K  and  in  Figure  2  for  the 
nonequilibnum  case  with  translational  temperature,  T^  =  14,000  K,  rotational  temperature  T 
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=  5000  K,  and  vibrational  temperature,  =  1000  K.  The  QCT  results  are  also  shown  in  the 
figures.  As  expected,  the  agreement  is  quite  good.  Thus,  the  increased  production  of  NO  in  the 
bow  shock  case  that  has  been  observed  with  the  implementation  of  the  GCE  model10  is  likely 
due  to  the  larger  overall  rate  constant  for  Reaction  (1)  rather  than  the  dependence  on  internal 
energy  state  distribution. 


Ev  (eV) 


Figure  1 

Vibrational  Distribution  of  Reacting  N2  Molecules.  =  Trot  =  5000  K 


Figure  2 

Vibrational  Distribution  of  Reacting  N2  Molecules. 
Ttnns  »  14000  K,  Trot  =  5000  K,  =  1000  K 
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DSMC  DISSOCIATION  MODELS 

Collision-induced  dissociation  reactions  of  diatomic  molecules  are  often  considered  to  be 
vibrationally  favored,  yet  unambiguous  information  on  the  degree  of  favoring  is  scant  The 
number  of  methods  that  have  been  used  in  the  context  of  DSMC  simulations  is  fairly  limited  In 
addition  to  the  TCE  and  EAE  models  described  above,  Haas  and  Boyd11*12  introduced  the 
vibrationally  favored  dissociation  (VFE>)  model  which  modifies  the  TCE  model  and  allows  for 
favoring  of  vibrational  energy  (Ev)  in  the  dissociation  probability  through  an  adjustable 
Paran^e.I\  Ivanov,  et  al.*  have  implemented  a  dissociation  model  in  their 
DSMC  code  that  uses  probabilities  based  on  state-specific  rate  constants  that  include  a  certain 
(fixed)  amount  of  vibrational  favoring. 

Macheret  and  Rich11-12  have  presented  an  analytical  threshold  line  model  which  includes 
ptr,  br,  and  Ev  as  parameters.  Recently,  Boyd13  has  adapted  this  model  for  a  DSMC 
importation,  using  one  probability  expression  for  dissociation  from  high  vibrational  levels 
^different  (reduced)  probability  expression  for  low  vibrational  levels.  The  reduction  in 
probability  for  lower  levels  means  that  this  model  is  vibrationally  favored.  For  the  case  of  N?-N? 
presented  in  Reference  13,  the  transition  occurs  around  v  =  6.  This  model  has  one  adjustable 
parameter  that  is  used  to  match  the  equilibrium  dissociation  rate  constant  at  some  temperature 
but  uses  no  adjustable  parameter  relating  to  the  degree  of  vibrational  favoring.  Boyd  has  shown 
mat  the  prediction  of  the  threshold  line  model  for  the  equilibrium  dissociation  rate  constant  as  a 
function  of  temperature  is  reasonable.  The  threshold  line  model  as  described  in  Reference  13  has 
been  implemented  in  the  code,  but  the  effect  of  mathematical  singularities11’12  near  certain 
values  of  vibrational  energy  has  been  observed  to  be  significant.  Use  of  the  foil  threshold 
model  >  rather  than  the  simplified  model13  reduces  the  effect  of  the  singularities  and  limits  the 
spurious  peak  that  occurred  in  the  region  of  the  transition  vibrational  level  as  can  be  seen  in 
Figure  3  for  the  case  of  H2  dissociation.  Details  of  the  implementation  are  given  in  Reference  14. 

lhe  H2  threshold  line  model  results  shown  in  Figures  4  and  5  here  all  used  the  corrected  version 
or  the  model. 


♦Ivanov,  M.S.,  Gimelshein,  S.F.,  Kashovsky,  A.V.,  and  Markelov,  G.N.,  private  communication. 
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Figure  4 

Comparison  of  QCT  Results  for  H2  +  Ar  with  DSMC  Dissociation  Models 

Figure  5 


Ev  (cm  “1  ) 

Threshold  Dissociation  Model  for  H2  +  H;  Comparison  with  QCT  Results 

The  works  of  Levine,  et  alA$  describe  in  detail  a  compact  method  to  quantify  the 
enhancement  of  chemical  reactions  due  to  reactant  internal  energy  using  a  single  parameter,  X. 
When  X  =  0,  then  there  is  no  vibrational  favoring,  and  the  reaction  probability  is  increased  by  all 
types  of  energy  equally.  The  magnitude  of  X  indicates  the  degree  of  vibrational  favoring  An 
approach  to  modeling  vibrationally  favored  dissociation  based  on  the  formalism  of  Levine  has 
been  used  in  several  studies  15'| 8  for  continuum  simulations  of  nonequilibrium  gases.  For  particle 
simulations,  Marriott,19  Gallis,20  and  Koura21  have  presented  adaptations  of  the  Levine 
approach;  all  of  these  include  an  adjustable  parameter,  X,  that  indicates  the  degree  of  vibrational 
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favoring,  but  formulated  in  different  ways.  Koura  bases  his  dissociation  model  on  the  weak 
vibrational  bias  (WVB)  equation  of  Reference  20  and  recommends  X  =  2  for  02-argon. 

The  method  of  displaying  the  vibrational  distribution  of  the  reacting  molecules  shown  in 
Figure  3  (and  in  Figures  1  and  2  for  the  case  of  an  exchange  reaction)  is  useful  in  visualizing  the 
differences  in  dissociation  models.  This  has  been  used  by  Boyd,13  who  compares  the  probability 
distribution  for  N2-N2  dissociation  in  equilibrium  at  10,000  K  for  the  VFD  (using  an  estimated 
value  of  <j>)  and  threshold  line  models.  The  results  demonstrate  that,  even  though  two  models  may 
reproduce  the  same  equilibrium  dissociation  rate  constant  and  both  are  vibrationally  favored, 
they  may  have  very  different  behaviors  at  the  state-specific  level.  Unfortunately,  state-specific 
information  is  not  available  for  N2  dissociation  (or  for  O2),  so  choosing  the  most  realistic  model 
is  difficult.  A  previous  study6  addressed  some  of  the  difficulties  in  using  shock  tube  data  to 
distinguish  between  chemistry  models  with  different  amounts  of  vibrational  favoring.  Shock  tube 
measurements  are  typically  of  macroscopic  parameters  (such  as  a  nonequilibrium  dissociation 
rate  constant  and  dissociation  induction  time),  which  are  not  very  sensitive  to  the  probability 
distribution  and  do  not  provide  the  state-specific  information  desired  to  validate  a  DSMC  model. 

To  begin  to  address  these  issues,  this  study  made  use  of  the  state-specific  dissociation 
results  that  are  available  from  QCT  calculations  for  H2.  Truhlar,  et  al.22  have  obtained 
vibrational  state  (rotationally  averaged)  dissociation  rate  constants  for  H2  dissociation  in 
collisions  with  argon  at  4500  K.  Also,  Martin  and  Mandy23  obtained  dissociation  rate  constants 
for  every  vibration-rotation  level  of  H2  colliding  with  H  for  the  entire  temperature  range  of  450- 
45,000  K.  In  order  to  obtain  rotationally  averaged  results  for  the  H2-H  case,  the  rate  constants  for 
each  rotational  level  within  a  given  vibrational  level  have  been  appropriately  weighted 
(according  to  degeneracy  and  Boltzmann  fraction)  and  summed  over. 

For  the  first  comparison,  an  equilibrium  case  for  H2-argon  at  4500  K  was  calculated  using 
three  DSMC  dissociation  models:  EAE,  WVB  as  described  in  Reference  20,  and  threshold  line 
model  Figure  4  compares  the  QCT  results22  with  DSMC  calculations  for  the  EAE,  threshold, 
and  WVB  models'  predictions  of  the  vibrational  distribution  of  dissociating  H2.  The  WVB  model 
result  is  shown  for  X  —  2.  The  model  was  also  run  for  X  =  A,  which  is  not  shown.  Increasing  X 
causes  the  distribution  to  be  peaked  more  toward  higher  v  and  the  two  values  bracket  the  QCT 
points  quite  well.  (Reference  22  states  that  a  fit  of  their  results  to  the  WVB  equation  yields  a 
value  ofX  =  2.2,  but  the  fit  in  this  study  using  Truhlar's  result  yields  X  =  2.8).  The  probability 
distribution  shape  of  the  EAE,  with  no  vibrational  favoring,  is  significantly  different  than  the 
QCT  result.  The  WVB  model  reproduces  the  correct  shape  when  one  uses  the  proper  value  of  X, 
but  there  is  the  problem  of  how  one  is  to  obtain  X  for  any  molecule.  The  threshold  line  model, 
with  no  adjustable  parameter,  shows  very  favorable  agreement  with  the  QCT  result. 

Very  similar  results  are  shown  for  H2-H  dissociation  in  Figure  5  for  the  threshold  model 
compared  with  QCT  results23  at  4500  K.  The  EAE  and  WVB  results  for  this  case  are  not  shown, 
but  give  basically  the  same  shape  for  H2-H  as  for  H2-Ar.  The  QCT  results  for  H2-H  provide 
some  information  on  how  the  probability  distribution  should  change  with  temperature.  The  state- 
specific  QCT  results21  have  been  converted  to  thermal  rotationally  averaged  values  for  three 
temperatures:  1000  K,  4500  K,  and  20,000  K.  The  vibrational  probability  distributions  for  the 
three  temperatures  are  compared  in  Figure  6.  As  can  be  seen,  the  effect  of  vibrational  favoring 
becomes  less  pronounced  as  the  temperature  increases.  This  trend  is  predicted  in  the  works  of 
Levine4’5  who  formulate  an  expression  for  vibrational  favoring  as  exp(-A,*Ev/kT).  This 
expression  has  been  shown  to  be  realistic  for  exchange  reactions  where  a  translation-rotation 
temperature  exists,  but  is  not  suited  to  DSMC  calculations.  The  WVB  model  does  not  include 
any  provision  for  a  reduction  in  vibrational  favoring  at  higher  temperatures,  and  so  may  be 
expected  to  be  realistic  over  some  range  of  temperatures,  but  less  so  as  the  range  becomes  very 
wide.  Preliminary  results  using  the  threshold  line  model  indicate  this  model  may  produce  a 
reduction  in  vibrational  favoring  as  temperature  is  increased,  but  further  studies  are  needed. 
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CONCLUSIONS 

Comparisons  with  QCT  results  for  the  first  Zel’dovich  reaction  confirm  that  the  EAE 
chemistry  model  gives  a  very  realistic  dependence  on  internal  energy  states  for  a  reaction  that  is 
not  vibrationally  favored.  Naturally,  it  is  not  to  be  expected  that  the  probability  distribution 
predicted  by  the  EAE  model  will  be  accurate  in  the  case  of  an  exchange  reaction  that  is  known  to 
be  strongly  vibrationally  favored,4’5  such  as  HC1  + 1. 

For  the  case  of  H2  dissociation,  the  WVB  model  with  X  =  2  to  4  and  the  threshold  line 
model  give  qualitative  agreement  with  QCT  results,  showing  the  correct  type  of  vibrational 
favoring  effect.  It  should  be  noted,  however,  that  the  H2  molecule  is  not  necessarily 
representative  of  other  diatomics.  Its  widely  spaced  rotational  and  vibrational  levels  cause 
quantum  effects  to  play  a  more  important  role.  Accurate  state-specific  dissociation  probabilities 
or  rate  constants  are  needed  for  other  diatomic  molecules  to  further  assess  DSMC  models. 
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PART  III:  DSMC  COMPARISONS  TO  ROTATIONAL  TEMPERATURE 
MEASUREMENTS  IN  JET  EXPANSIONS  WITH  FINITE  BACKGROUND 
PRESSURES 


INTRODUCTION 

A  comparison  has  been  made  between  Direct  Simulation  Monte  Carlo  (DSMC) 
calculations  and  experimental  measurements  of  rotational  temperature  in  a  nitrogen  jet 
expanding  into  a  finite  background  pressure  in  which  a  distinct,  but  non-continuum,  Mach  disk  is 
formed.1  In  the  work  presented  in  this  report,  corrections  and  improvements  to  various  parts  of 
the  simulation  have  been  made,  a  new  stagnation  region/orifice  calculation  has  been  carried  out 
to  supply  a  new  startline  to  the  expansion  flow  calculation,  and  an  investigation  of  the  effects  of 
a  revised  rotational  collision  model  at  low  temperatures  has  been  accomplished. 

The  experimental  results  that  are  used  for  comparison  to  the  simulations  were  obtained  at 
the  University  of  California  at  Berkeley  Rarefied  Gas  Wind  Tunnel.  In  those  experiments, 
electron  beam  fluorescence  from  the  (0,0)  vibrational  band  of  the  first  negative  system  of  the 
nitrogen  ion  was  used  to  spectroscopically  probe  the  centerline  of  a  heated  nitrogen  jet  expanding 
into  a  chamber  maintained  at  a  low  but  finite  background  pressure.  Centerline  rotational 
temperatures  were  determined  for  a  number  of  conditions.  A  complete  description  of  the 
experimental  procedures  and  results  is  presented  by  Gochberg2  and  a  summary  is  given  in 
Reference  1. 


ROTATIONAL  RELAXATION  MODELS 

The  prediction  of  the  rotational  temperature  of  a  gas  in  an  expansion  flow  or  compression 
shock,  in  which  large  deviations  from  equilibrium  with  the  translational  mode  of  the  gas  can 
occur*  particularly  sensitive  to  the  treatment  of  collisional  energy  exchange  in  the  simulation. 
In  Reference  1,  it  was  demonstrated  that  a  temperature  dependent  rotational  collision  number 
model  produced  significantly  different  rotational  temperature  profiles  in  a  jet  expansion  than  that 
of  a  constant  collision  number  model.  In  addition,  the  viscosity  coefficient  used  in  the  variable 
hard  sphere  collision  model  was  also  shown  to  affect  the  calculated  rotational  temperatures 
Gochberg  and  Haas  •  have  examined  the  effects  of  rotational  relaxation  rate  models  on  the 
rotational  temperature  m  vacuum  free  jets  using  an  analytical  method  that  involved  numerical 
integration  of  the  governing  rotational  relaxation  rate  equations.  In  comparing  the  theoretical 
results  to  experimental  measurements,  they  were  unable  to  recommend  an  optimum  rotational 
relaxation  model,  or  combination  of  models,  that  would  give  a  good  match  to  all  of  the 
experimental  data,  due  primarily  to  the  large  degree  of  scatter  in  the  data.  They  did  conclude  that 
the  Parker  model  does  not  accurately  predict  the  rotational  temperature  at  low  translational 
temperatures,  due  to  the  fact  that  this  model  predicts  rotational  collision  relaxation  times  that  go 
to  zero  with  decreasing  temperatures.  Quantum  effects  are  known  to  become  significant  in 
rotational  relaxation  at  low  temperatures,  and  the  model  of  Lebed’  and  Riabov6  is  one  example  of 
.  of  accounting  for  these  effects.  For  temperatures  below  approximately  125  K,  quantum 
effects  increase  the  rotational  relaxation  time  by  an  order  of  magnitude  or  more  compared  to  the 

Parker  model,  slowing  the  equilibration  of  the  rotational  mode  to  the  translational  mode  of  the 
gas. 


DSMC  SIMULATION 

The  key  features  that  must  be  addressed  in  a  DSMC  simulation  of  this  problem  are  the 
plenum,  jet,  and  vacuum  chamber  boundary  conditions,  and  the  rotational  energy  relaxation 
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model. 


For  the  relatively  high  plenum  pressures  of  the  Gochberg  experiments  it  is  not  possible  to 
do  a  complete  DSMC  calculation  for  the  entire  plenum-jet  system.  Therefore,  in  this  study,  the 
calculations  were  broken  up  by  starting  the  expansion  jet  DSMC  simulation  from  a  position  just 
upstream  of  the  orifice  exit  plane.  Providing  an  inflow  startline  at  this  position  is  problematic. 
The  orifice  used  in  the  experiments  had  an  1/d  of  approximately  0.6,  and  so  was  not  an  ideal 
sonic  orifice  with  Mach  1  at  the  exit  plane.  Analysis  of  these  types  of  tube  flows  have  shown 
that,  dong  the  centerline  of  die  flow,  sonic  conditions  will  occur  slightly  upstream  of  the  exit 
plane.  ’  An  isentropic  expansion  from  the  measured  experimental  stagnation  conditions  to  sonic 
conditions  (Mach  1)  at  a  position  0.1  mm  (0.12  diameters)  upsteam  of  the  exit  plane  was  used.  A 
“plug”  flow  (no  boundary  layer)  was  assumed  at  this  position  to  define  the  inflow  startline  for 
the  jet  expansion  flow  calculation.  Simulations  to  determine  the  effects  on  centerline  jet 
properties  of  a  full  boundary  layer  were  carried  out,  and  the  results  indicated  a  finite  but  small 
influence  of  the  boundary  layer  structure  as  long  as  the  centerline  inflow  startline  parameters 
were  identical. 

It  would  be  advantageous  to  use  a  full  calculation  of  the  flow  to  better  define  the  conditions 
at  the  inflow  startline  position.  DSMC  simulation  of  the  flow  into  the  orifice  from  the  stagnation 
region  at  the  relative  high  pressures  of  the  experiment  (250-550  torr)  is  difficult  even  on  the 
largest  computers.  Navier-Stokes  calculations  are  also  difficult  due  to  the  fast  expansion 
characteristics  of  this  flow,  as  well  as  the  difficulty  of  defining  the  boundary  conditions.  The 
stagnation  region  DSMC  calculation  used  in  Reference  1  had  extremely  large  cells  compared  to 
the  , local  mean  free  path  of  the  gas,  and  so  was  unreliable.  The  pressure-matching  method  of 
Serikov  was  implemented  and  used  for  the  inflow  boundary  in  the  stagnation  region.  This 
allowed  a  calculation  to  be  made  using  a  computational  region  only  1  orifice  diameter  upstream 
of  the  entrance  to  the  orifice.  A  calculation  of  the  flow  from  the  stagnation  region,  through  the 
orifice,  and  out  a  small  distance  into  the  jet  expansion  was  accomplished  on  a  large  parallel 
computer,  allowing  cell  sizes  to  be  comparable  to  the  local  mean  free  path.  This  calculation 
produced  sonic  conditions  well  upstream  of  the  exit  plane  and  a  much  higher  translational 
temperature  and  velocity  than  the  isentropic  expansion  prediction.  This  might  be  due  to  the 
effects  of  the  hot  wall  on  the  flow  (all  walls  in  the  calculation  are  set  to  the  stagnation  region 
temperature  of  753  K).  Table  1  presents  the  centerline  values  of  the  flow  parameters  for  die  sonic 
and  DSMC  startlines  at  a  position  0. 1  mm  upstream  of  the  orifice  exit  plane. 


Type 

Jiautc  jl*  lunui 

Mach 

Number 

w  aiarume  r  tow  ran 
Number  Density 
(m  ) 

ameters  on  uente 
Temperature 
(K) 

mine 

Velocity 

(m/s) 

DSMC 

Sonic 

1.15 

1.0 

2.49x10“ 
2.858  xlO24 

595 

477 

580 

445 

As  in  Reference  1,  collision  energy-dependent  rotational  collision  numbers  based  on  the 
Parker  model,  as  generalized  for  variable  hard  sphere  (VHS)  molecules  by  Boyd,10  was  used 
(with  Chung’s  correction  factor11),  as  well  as  the  appropriate  low  temperature  value  of  w.  The 
corresponding  DSMC  probabilities  of  relaxation  were  based  on  the  methods  of  Haas,  et  al.n  This 
method  produces  N2  -  N2  collisional  rotational  relaxation  rates  in  the  DSMC  simulation  about  a 
factor  of  2  higher  than  the  P  =  1/Z  assumption. 

An  error  in  the  treatment  of  the  chamber  background  gas  boundary  condition  in  the  DSMC 
calculations  was  corrected,  which  eliminated  the  large  difference  in  the  position  of  the  rotational 
temperature  shock  between  the  DSMC  prediction  and  the  experimental  data  that  was  reported  in 
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Reference  1 .  Also,  various  coding  errors  in  the  handling  of  internal  energy  transfer  in  collisions 
were  corrected.  This  prevented  a  steady  loss  of  rotational  energy  that  was  previously  occurring. 

RESULTS 

All  of  the  calculations  presented  here  are  for  the  conditions  of  Case  1  in  Gochberg’s2 
experiments:  0.8125  mm  diameter  orifice,  T0  =  753  K,  P0  =  351.6  torr,  Toe  =  297.8  K,  Poo  =  48 
mtorr.  The  centerline  rotational  temperature  profiles  for  the  isentropic  “plug”  flow  and  DSMC 
calculation  startlines  are  shown  in  Figure  1 .  The  isentropic  startline  gives  a  closer  match  to  the 
data,  but  both  calculations  produce  rotational  temperatures  higher  than  the  data,  even  in  the  post- 
shock  region.  It  should  be  pointed  out  that  the  gas  rotational  temperature  is  completely 
equilibrated  to  the  translational  temperature  in  the  post  shock  region  for  all  the  simulations.  Also 
shown  in  the  figure  is  the  result  for  a  constant  rotational  collision  number  model  using  a  value  of 
.  that  parameter.  As  was  found  in  Reference  1,  the  constant  collision  number  model  predicts 
significantly  different  rotational  temperatures  than  the  variable  rotational  collision  model,  in  a 
direction  farther  from  the  experimental  data.  Spatially  extended  calculations  show  that  the  post 
shock  conditions  slowly  approach  the  chamber  conditions  over  a  distance  of  hundreds  of  orifice 
diameters  downstream  from  the  shock  position.  Because  the  DSMC  calculation  startline  may  still 
be  problematic,  the  isentropic  startline  was  used  in  all  of  the  analyses  reported  here. 


Different  schemes  were  used  to  define  the  outflow  boundary  in  the  DSMC  calculation.  The 
boundary  of  the  calculation  was  also  placed  at  different  distances  downstream.  It  was  found  that 
as  long  as  the  calculational  region  extended  to  about  twice  the  distance  to  the  shock,  the 
boundary  treatment  and  calculational  boundary  position  did  not  affect  the  quantitative  structure 
of  the  shock. 

A  portion  of  the  expansion  flow  upstream  of  the  shock  attains  translational  temperatures 
beiow  the  125  K  level  at  which  the  Parker  rotational  relaxation  model  is  known  to  be  incorrect. 
•  i-'r-  j  ^  implementation  of  the  Lebed’  and  Riabov6  rotational  relaxation  model  a 

scheme  was  used  whereby  the  rotational  relaxation  collision  number  was  set  to  a  value 
of  10  tor  collisions  m  which  the  collisional  energy  was  below  an  equivalent  temperature  of  125 
K.  from  these  results,  it  was  demonstrated  that  the  influence  of  the  low  temperature  region 


31 


rotational  relaxation  process  on  the  rotational  temperatures  in  the  shock  region  was  negligible.  A 
correction  to  the  rotational  relaxation  model  at  low  temperatures  is  therefore  not  necessary  in 
comparing  the  predicted  and  measured  rotational  temperatures  in  the  shock  region. 

To  investigate  the  effects  of  changes  in  the  calculational  boundary  conditions  on  the  shock 
structure,  a  series  of  computations  were  made  for  a  range  of  pressure  ratios  by  varying  the 
number  density  at  the  inflow  startline  and  at  the  outflow  vacuum  chamber  boundaries.  The 
values  of  the  parameters  are  shown  in  Table  2  and  the  results  are  shown  in  Figure  2.  An 
increased  pressure  ratio  (Po/P  )  causes  the  shock  position  to  move  downstream  as  expected.  The 
post  shock  rotational  temperature  is  not  affected  significantly  by  the  pressure  ratio  over  the  range 
investigated.  The  peak  height  of  the  rotational  shock  is  not  affected  by  the  startline  number 
density,  but  is  dependent  on  the  chamber  density.  The  slope  of  the  shock  is  also  affected 
significantly  by  the  chamber  density.  The  significant  change  in  the  shock  position  with  a  small 
(10%)  change  in  the  startline  number  density  points  out  the  importance  of  an  accurate  startline 
for  these  conditions. 


Table  2.  Chamber  and  Startime  Number  Density  Values 


Chamber  Number 

Inflow  Startline 

Density  (m'3) 

Number  Density  (m3) 

Baseline 

1.556x10"* 

2.858  x  10"4 

Low  Chamber 

1.0  xlO21 

2.858  x  1024 

High  Startline 

1.556  xlO21 

3.1438  x  1024 

Figure  2 

Effects  of  Pressure  Ratio 


To  further  understand  the  effects  of  the  inflow  startline  on  the  predicted  rotational 
temperature  shock  structure,  two  runs  were  made  using  inflow  startline  velocity  and  temperature 
10%  higher  than  the  baseline  values  for  the  sonic  “plug”  startline.  The  results  are  shown  in 
Figure  3.  These  small  changes  in  inflow  startline  velocity  and  temperature  produced  a  significant 
increase  in  the  peak  and  post-shock  rotational  temperature  and  tended  to  flatten  the  profile.  The 
rotational  temperature  shock  position  was  moved  downstream  for  the  higher  velocity  case,  but 
did  not  change  for  the  higher  temperature  condition.  Once  again  it  is  seen  that  small  changes  in 
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the  inflow  startline  can  produce  significant  changes  in  the  predicted  shock  structure. 


Figure  3 

Effects  of  Startline  Velocity  and  Temperature 


CONCLUSIONS 


The  improvements  to  the  DSMC  rotational  relaxation  collision  model  implementation  that 

nrl-fii  been  ma^e  smce  1116  study  m  Reference  1  was  completed  have  changed  the  predicted 
profiles  of  rotational  temperatures  m  the  expansion  jet  shock.  When  compared  to  the 

da1f /or  Gochberg  s  Case  1,  the  improvement  is  mixed.  The  position  of  the  shock 
more  closely  matches  the  experimental  results,  but  the  peak  height  and  post-shock  temperature  is 
now  significantly  greater  than  the  data  values.  F  temperature  is 


„  F°r  the  low  temperatures  present  in  the  initial  part  of  the  jet  expansion,  it  is  possible  that  the 
Sff  hard  sphere  (GHS)  collision  model  (Hassan  and  Hash13),  which  includes  an  attractive 
component  of  the  molecule-  potential,  may  be  more  accurate  than  the  repulsive  power-law 
present  m  the  VHS  model.  For  example,  in  an  equilibrium  gas  at  T  =  50  K,  the  GHS  collision 

i  awr??/ 1S  6i s  VHS  ^ue  («$«  GHS  311(1  VHS  parameters  given  in  Reference 
).  This  would  tend  to  increase  the  eqmlibration  of  Tr  to  Tt  in  the  expansion  region.  The  GHS 

Chnm,^  unplemented  usmg  the  rotational  relaxation  methods  discussed  in  Hash,  et  a/.,14  and 
Because  of  a  severe  increase  m  computational  time  requirements  for  the  collision 
selection  procedures,  it  was  not  been  possible  to  obtain  solutions  for  this  flow  with  the  GHS 
model.  As  an  alternative,  the  variable  soft  sphere  (VSS)  model  (Koura  and  Matsumoto16)  with 
the  appropriate  low  temperature  parameters  has  been  considered  but  not  yet  implemented. 
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APPENDIX 


DSMC  Collision  Procedure 

The  basic  Direct  Simulation  Monte  Carlo  (DSMC)  methodology  is  well  described  in 
Bird.-*-  The  present  study  focused  on  the  collision  procedure,  specifically  the  degree  to  which 
it  can  model  certain  types  of  nonequilibrium  physics.  The  collision  processes  considered  are 
summarized  and  consistent  implementation  of  models  into  the  DSMC  algorithm  is  discussed, 
including  key  differences  between  the  present  approach  and  those  in  the  literature. 

Molecular  Model 


The  present  study  was  limited  to  atomic  and  diatomic  species  in  the  ground  electronic 
state.  For  diatomic  species,  the  internal  modes  (rotation  and  vibration)  are  nominally  con¬ 
sidered  uncoupled,  though  some  portions  of  this  analysis  could  be  applied  to  coupled  states 
(including,  e.g.,  quasi-bound  rotational  levels). 

Rotational  Mode.  The  rotational  mode  is  described  by  a  bounded  version  of  the 
quantized  rigid  rotator  model  proposed  by  Boyd.^  The  rotational  energy  at  level  j,  0  <  j  < 
jo,  given  by  er(j)  =  j(j  +  1  )k6r,  where  9r  is  the  characteristic  rotational  temperature  (see, 
e.g.,  Bird,*-  Appendix  A)  and  k  is  the  Boltzmann  constant.  The  maximum  bound  level  is 
given  by 

h  =  {-1  +  [1  +  (4D/Wr)]>'2}/2  (1) 

where  D  is  the  dissociation  energy  (measured  from  the  ground  vibrational  state). 

The  normalized  equilibrium  distribution  of  eT  at  temperature  T  is  defined  as 


I'b  ti-.T) 

=  fs(r,T)/Qr(T ) 

(2) 

feti-.T) 

=  (2j-  +  1)  exp  (-f.JkT) 

j=3D 

(3) 

Qr(T ) 

=  E/bU;T) 

j=o 

(4) 

Replacing  the  sum  in  the  partition  function  Qr  by  integration  gives 


S'bU;T)~ 


0r,„,  ,  ,,  exp  (-tr/kT) 
T(2j  +  1  1  -  exp  (—D/kT) 


(5) 


This  reduces  to  Equation  (2)  of  Boyd  for  the  case  kT  <C  D.  Initial  rotational  energies  for 
particles  are  defined  by  sampling  from  this  distribution. 


The  code  outputs  a  rotational  temperature  calculated  from  the  sampled  average  energy 
er  by  assuming  the  equilibrium  distribution  of  Equation  (3)  and  solving  for  the  parameter 
Tr  in 


-  _  go DjMjjTr) 
€r  Qr(Tr) 


(6) 
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It  should  be  noted  that  the  rotational  distribution  function  represented  by  diatoms  in  the 
simulation  need  not  be  equilibrium,  and  this  temperature  should  be  interpreted  as  simply 
an  estimate  of  the  average  energy  present  in  the  rotational  mode. 

Vibrational  Mode.  The  vibrational  mode  is  described  by  the  bound  quantized  an- 
harmonic  oscillator  model  of  Bergemann  and  Boyd.^  Bound  vibrational  energy  levels  v,  0  < 
v  <  vD,  of  energy  ev(v)  are  nominally  generated  from  the  Dunham  expansion  using  the 
spectroscopic  constants  available  in  Huber  k  Herzberg.4  (Note:  the  zero  point  energy  was 
subtracted  to  obtain  e„(0)  =  0). 

The  normalized  equilibrium  distribution  of  ev  at  temperature  T  is  defined  as 


f’s^T) 

=  /s(< r,T)/Q,(T) 

(7) 

fs(v;T) 

=  exp  (-c„/kT) 

V—VD 

(8) 

Q,(T) 

=  Y.  fe(v;T) 

v—0 

0) 

Initial  vibrational  energies  for  particles  are  defined  by  sampling  from  this  distribution.  Note 
that  for  anharmonic  levels,  the  inversion  of  the  distribution  function  to  assign  particle  vibra¬ 
tional  energy  as  used  by  Bird1  (Appendix  C,  and  implemented  in  the  sample  programs)  is 
not  valid  and  will  result  in  incorrect  distributions  (with  the  error  increasing  with  T).  That 
method  is  valid,  however,  for  harmonic  oscillators  (constant  energy  gap  ev(v)  =  vkdv). 

The  code  outputs  a  vibrational  temperature  Tv  calculated  from  the  sampled  energy  ev 
analogously  to  that  for  Tr. 

Collision  Procedure 


The  collision  procedure  is  based  on  the  serial  redistribution  method  of  Bird1  (Chap.  5.5), 
as  extended  by  Haas,  et  al.^  A  schematic  of  the  basic  procedure  is  given  in  Figure  Al.  The 
number  of  pair  selections  per  time  step  under  the  No  Time  Counter  (NTC)  model  (Bird,1 
Chap.  11.2),  is  based  on  the  maximum  value  of  the  product  of  the  elastic  cross  section  aei(et ) 
(here  given  by  the  Variable  Hard  Sphere  (VHS)  model),  and  the  relative  speed  g  defined  by 
the  relative  translational  energy  et  =  mrg2/ 2,  mT  being  the  reduced  mass  of  the  collision  pair, 
and  will  thus  recover  the  elastic  relaxation  rate.  Collision  pairs  are  randomly  selected  and 
are  then  accepted  with  the  probability  P  =  crelg/(aeig)max  by  comparing  P  with  a  uniform 


random  number  R,  0  <  R  <  1.  The  code  allows  for  the  species  to  be  combined  into  ‘groups’ 
(Bird,1  Chap  11.2),  with  a  (aeig)max  stored  for  each,  to  improve  acceptance  efficiency  for 


cases  where  large  disparities  between  nominal  aei  or  g  arise  between  species. 


Details  of  Collision  Events.  All  accepted  collision  pairs  undergo  redistribution  of 
translational  energy  as  the  last  process  of  the  collision  procedure.  The  available  inelastic 
or  chemical  reaction  collision  events  are  dependent  on  the  type  of  the  species  pair  selected. 
For  atom-atom  collisions,  only  a  recombination  reaction  is  potentially  available.  For  diatom- 
atom  collisions,  a  dissociation  and,  possibly,  an  exchange  reaction,  along  with  vibrational- 
translational  (V-T)  and  rotational-translational  (R-T)  relaxation,  are  available.  Diatom- 
diatom  collisions  are  similar  to  those  for  diatom-atom  pairs;  however,  the  exchange  reaction 
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is  not  available.  Figure  A2  shows  the  process  by  which  possible  collision  outcomes  are 
modeled  for  the  case  of  a  diatom-diatom  pair,  denoted  A2 — B2,  which  has  been  accepted  for 
collision  via  the  logic  of  Figure  Al.  The  figure  shows  only  the  case  where  the  first  particle 
chosen  for  testing  is  A2.  The  case  for  B2  chosen  first  is  completely  analogous.  Note  that  the 
diatoms  need  not  be  homomolecular. 


The  serial  procedure  includes  the  ‘particle  selection  prohibiting  multiple  relaxation’  logic 
of  Haas,  et  al.  to  calculate  the  probability  of  each  event.  This  method  is  an  attempt  to 
improve  the  numerical  treatment  of,  rather  than  the  physical  modeling  of,  several  dependent 
collision  processes.  The  logic  analytically  relates  the  molecular  probability  Px  of  a  relaxation 
event  of  mode  x  and  the  macroscopic  average  value  (typically  described  by  the  input  collision 
number  Zx)  under  equilibrium  conditions.  The  collision  number  is  defined  as  ratio  of  the 
mode  x  relaxation  time  to  an  elastic  collision  time,  Zx  ~  rx/rc,  so  approximately  Px  ~  1/ZX. 
The  Bird  serial  NTC  methodology  determines  the  total  number  of  collisions  from  the  elastic 
cross  section.  The  Haas,  et  al.5  logic  increases  the  probabilities  for  inelastic  events  (P  >  l/Z) 
to  obtain  the  proper  relaxation  rate  (characterized  by  Z')  of  each  mode.  Using  this  approach, 
the  event  cross  section  <jx  ~  Pxcei  can  be  identified.  Note  the  requirement  that  a  microscopic 
probability  P  or  cross  section  <rx  be  related  to  (ill-defined)  macroscopic  quantities  such  as  Z 
is  a  common  weakness  of  the  majority  of  DSMC  collision  models. 

In  the  present  study,  the  Haas,  et  al.5  logic  has  been  applied  only  to  the  internal  relaxation 
(vibration  followed  by  rotation)  portion  of  the  complete  serial  procedure  of  Figure  A2.  The 
formulae  for  the  P  are  dependent  on  the  ordering  of  modes  in  the  serial  procedure  (see,  e.g., 
Appendix  B  of  Haas,  et  al.).  For  the  case  of  vibrational  relaxation  in  the  diatom-diatom 
collision  of  Figure  A2,  each  probability  is  dependent  on  the  relaxation  rate  values  of  both  of 
the  colliding  species,  and  the  PV2  of  the  second  diatom  (B2  in  Figure  A2)  is  dependent  on 
the  value  calculated  for  the  first  diatom  A2  tested  (but  not  accepted)  PVl : 

pv  1  =  root  of  [J*  +  (X2  -  Xi~  2)Pn  +  2A:  =  0]  (10) 

P  -  *2  ,  , 

(l-P.,/2)  <n> 


where 


X\  = 


Ct  Cu  1 


(12) 


C t  4, 

with  all  quantities  in  Xx  those  for  impact  of  species  2  on  target  diatom  1,  and  similarly 
for  X2.  The  values  of  Z  are  as  discussed  below.  The  £  terms  describe  the  degrees  of 
freedom  of  the  energy  modes  of  the  diatoms  appearing  in  the  Borgnakke-Larsen  procedure. 
The  translational  collision  energy  mode  (biased  by  the  NTC  collision  selection)  term  is 
Ct  =  5  -  2ub  and  the  vibrational  energy  mode  term  £»  =  C v(Tv)  =  2 ev/(kTv)  (Equation 
(11.28)  of  Bird),  based  on  the  sampled  cell  vibrational  temperature  of  the  target  species. 
As  shown  in  Figure  A2,  in  the  serial  (or  particle)  approach  prohibiting  multiple  relaxation, 
only  the  first  accepted  relaxation  event  is  carried  out.  For  example,  if  R  <  PVl)  then 


the  post-collision  vibrational  level  v1  of  A 2  is  calculated,  the  remaining  energy  e, 


e;  is 
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redistributed  into  the  translational  energy  (i.e.,  velocities)  of  the  colliding  pair,  and  the 
collision  is  concluded.  Otherwise  the  next  test  ( R  <  PV2)  is  made. 

Since  the  tests  for  rotational  relaxation  follow  those  for  vibrational  relaxation,  the  rota¬ 
tional  relaxation  probabilities  Pn ,  Pr2  become  dependent  on  those  for  vibration 

Pn  =  root  of  [P*  -  (2  -  Y2/Pv  +  Y1/Pv)Pn  +  2Y1/PV  =  0]  (13) 

D  _  y2 

T2  (1  —  Pn/2)PV  (14) 

Pv  =  (1  -  PV1)(1  -  PV2)  (15) 

where,  analogous  to  the  vibration  case 


(  Ct  +  Cr  1 

V  Ct  zr)l 


and  the  rotational  mode  contributes  (r  =  2. 


These  formulae  place  additional  mathematical  bounds  (which  may  or  may  not  have  a 
physical  basis)  on  the  allowable  Z.  In  contrast  to  the  discussion  in  the  paper  of  Haas,  et 
al., 5  the  method  has  been  found  here  to  be  valid  only  for  the  case  of  near-constant  Z  and 
will  not  accomodate  coupled  processes  for  which  each  possesses  a  Zx(ec).  Also,  note  that 
inclusion  of  additional  serial  collision  events  (e.g.,  as  shown  in  Figure  A2,  the  dissociation  test 
precedes  vibration)  means  that  the  coupling  between  processes  becomes  excessively  complex, 
and  this  method  does  not  appear  feasible  for  general  application.  The  null  collision  selection 
method  of  Koura6  is  less  mathematically  precise,  but  is  more  easily  generalizable  and  more 
physically  consistent.  It  appears  that  method  will  be  useful  in  future  studies,  but  will  not 
be  discussed  further  here. 


Cross  Sections  and  Probabilities.  Elastic  Cross  Sections.  Elastic  (and  thus  for 
the  nominal  Bird  method,  total)  cross  sections  are  based  on  VHS  values 


Oei{g)  =  aTe}etUB+l/2  (17) 

using  the  (constant)  reference  cross  section  ( aref )  and  exponents  in  the  viscosity  temper¬ 
ature  relationship  (/x(T)  ~  TUB)  tabulated  in,  e.g.,  Bird,'*'  Appendix  A.  Using  the  general 
Borgnakke-Larsen  procedure,  only  some  fraction  (~  \/Zx)  of  collision  pairs  accepted  for 
elastic  collisions  are  tested  for  chemical  reaction  or  inelastic  relaxation  events. 

Rotational  Relaxation.  Models  for  rotational  (specifically  rotational-translational  or 
R-T)  relaxation  are  generally  not  sophisticated  at  the  present  time.  All  available  DSMC 
models  apparently  are  based  on  the  theoretical  model  given  by  Equation  (45)  of  Parker,7 
which  predicts  qualitatively  the  experimentally  observed  temperature  dependence  of  the 
macroscopic  rotational  collision  number  Zt(T).  Boyd2  (see  also  Choquet8  and  Koura9)  have 
generalized  this  to  an  energy-dependent  form  Zr(ec)  for  use  in  DSMC  which  reproduces  the 
Parker  temperature  dependence  under  equilibrium. 
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Relatively  little  modeling  of  R-T  relaxation  was  done  using  the  present  DSMC  code  in 
this  part  of  the  effort.  For  cases  where  only  rotational  relaxation  has  been  allowed  in  the 
serial  procedure,  the  energy-dependent  form  ZT(ec)  given  by  Equation  (24)  of  Choquet® 
(rewritten  for  Zr )  has  been  used 


Zr  (ec) 


Zr°° 


1  +  4tt/3  {kT*/ecy/2  +  (tt2/4  +  2)(kT*/ec) 


(18) 


(note  the  slight  difference  in  the  final  constant  term  in  the  denominator,  chosen  to  match 
the  VHS  form  given  by  Equation  (18)  of  Hash,  et  al.^). 

For  cases  where  rotational  relaxation  was  preceded  by  the  test  for  vibrational  relaxation, 
the  temperature-dependent  model  Zr(T )  based  on  the  VHS  form  given  by  Equation  (18)  of 
Hash,  et  al.^  (see  also  Equation  (21)  of  Choquet^)  has  been  used 


Zr(T)  =  z, ; 


r(3  —  u>b) 


r(3  -  uB)  +  47r/3r(5/2  -  uB){T*/Ty /2  +  2(tt2/4  +  2)r(2  -  uB)(T*/T) 


(19) 


where  F()  is  the  gamma  function,  T  is  the  directly  sampled  translational  temperature  in  the 
cell  in  which  the  collision  partners  reside,  and  T*,  Z™  are  as  given  in,  e.g.,  Parker  7 

Note  that  using  the  serial  method  and  the  R-T  model  here,  only  a  single  particle  of 
the  selected  pair  is  allowed  to  relax  per  collision  event.  Consideration  of  the  possibility 
of  a  rotational-rotational  (R-R)  event  during  a  diatom-diatom  collision  would  require  the 
incorporation  of  relatively  complex  pairwise  probabilities  such  as  those  given  by  Equation 
(26)  of  Koura.^ 

Vibrational  Relaxation.  Several  models  of  vibrational-translational  (V-T)  relaxation, 
of  varying  levels  of  mathematical  and  physical  sophistication,  were  implemented.  Relatively 
comprehensive  comparisons  of  these  V-T  models  were  made  and  are  documented  as  part  of 
this  report. 


1.  ZV(T)  based  on  the  Millikan  &  White  (M&W)  correlation, ^  rv  =  exp  ( AT “V3  +  B)/(nkT), 
where  n  is  the  number  density,  T  is  the  cell  translational  temperature,  and  the  input 
constants  A  and  B  are  available  from  the  original  report  and  subsequent  compilations 
(see,  e.g.,  Moreau^). 

2.  The  M&W  value  including  the  empirical  Park^  correction  factor  which  limits  the 
minimum  allowable  rv  (and  thus  Zv),  ZV(T )  =  (rv  +  Tpark)/rc,  where  Tpark  is  some 
constant  C  or  C(T)  >  rc,  as  discussed  in  Park.  The  correction  factor  is  of  non- 
negligible  magnitude  only  for  relatively  high  temperatures  and  was  not  found  to  be 
important  in  the  present  work.  However,  by  fitting  results  to  the  more  sophisticated 
models  discussed  below,  improved  correction  factors  can  be  generated. 

3.  The  Bird^  (Chap  11.4)  form  of  the  M&W  model,  where  ZV(TC)  is  calculated  based  on 
a  temperature  corresponding  to  the  collision  energy,  Tc  =  ec/(Ct  +  Cv)k.  This  model, 
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using  the  M&W  rates  corresponding  to  Tc,  results  in  vibrational  collision  numbers  Zv 
one  or  more  orders  of  magnitude  too  small  (Choquet44)  and  thus  leads  to  excessive 
relaxation  rates.  This  error  can  be  seen  clearly  in  the  stagnation  streamline  flowfield 
of  Figure  12.47  of  Bird.  Those  calculations  use  the  serial  procedure,  similar  to  that 
shown  in  Figure  A2  here,  but  do  not  use  the  logic  of  Haas,  et  al.5  The  small  value  of 
Zv  means  the  majority  of  collisions  undergo  vibrational  relaxation,  and  relatively  few 
are  tested  for  rotational  relaxation.  This  gives  rise  to  the  non-physical  behavior  seen 
in  that  figure  that  Tv  >  Tr  through  the  shock  layer. 

4.  The  Abe4  ^  single-quantum-jump  harmonic  oscillator  model  valid  for  low  temperatures, 
calibrated  to  reproduce  the  M&W  relaxation  rate 

Pv,v—\  =  vP] io  (20) 

Pv,v+ 1  =  +  ,  (et>k6v)  (21) 

Pio  is  a  constant  probability  related  to  the  rate  constant  /c10  for  (v,  v')  =  (1, 0)  transi¬ 
tions  from  the  M&W  correlation  (see  also  Koura^) 

1  r°° 

*10  =  [1  -  exp  (— ev/T)\r.=l  P^f(9)d9  =  V(il2-^B)Pwl2  (22) 

The  net  probability  of  vibrational  relaxation  is  Pv  =  Pv>v-i  +PV,V+ 1.  Note  the  formulae 
given  by  Equation  (23)  and  (24)  of  Abe  assume  hard  sphere  collision  mechanics.  The 
formulae  here  are  generalized  to  the  VHS  model  (and  use  the  present  notation  coB  = 
w(Abe)  +  4/2>  an<4  t^ie  rate  consfant  ki0  =  fc10(Abe)/n)’  4r°r  a  ^ar<^  sphere,  Abe’s 
Equation  (24)  should  read  kito  =  nPiocrref[2kT/ his  x  is  extraneous. 


5.  The  energy-dependent  version  of  the  Abe  model  given  by  Equation  (22)  of  Abe  and 
the  related  Equations  (20)  and  (21)  of  Boyd,4^  e.g., 

(V)5+2w 

Pv,v+ 1  =  c(v  +  1)  ~~^2 —  exP  (-9*/ 9')  (23) 

where  C  and  g*  are  constants  and  the  post-collision  relative  speed  is  defined  by  e't  = 
mr(g')2/ 2  =  et  -  kdv. 

As  discussed  by  Choquet,44  these  formulations  lead  to  difficulties  with  equipartition 
of  energy  and  were  not  pursued  further. 

6.  The  Koura®  form  of  the  Heidrich,  et  al.4^  improvement  to  forced  oscillator,  impulsive 
transfer  semiclassical  (ITFITS)  multi-quanta-jump  harmonic  oscillator  model,  valid  for 
diatom-atom  collisions  over  a  wide  range  of  temperatures.  The  form  given  in  Koura 
apparently  uses  a  non-standard  polynomial  function,  and  for  clarity  Equation  (23)  of 
the  original  reference  is  used 


Pv,«  =  v!v'!((A F))v+v'  exp  (-(A E)) 


S5”')  ((A  E))-‘  \ 

h  «(»-!)!(»'-  OiJ 
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and  the  symmetrized  energy  gap  term  (A E))  is  as  defined  in  Equation  (22)  of  that 
reference.  The  total  probability  of  vibrational  relaxation  is  given  by  Pv  =  Pv,v> . 

Note  that  versions  of  the  ITFITS  model  applicable  to  anharmonic  oscillators10  and 
diatom-diatom  collision  partners1^  are  available,  along  with  a  similar  semiclassical 
model  for  V-V2(1  transfer,  but  have  not  yet  been  implemented. 


Dissociation.  Relatively  few  models  of  dissociation  are  available  which  incorporate  any 
degree  of  physical  realism.  It  is  expected  that  models  which  incorporate  an  (adjustable) 
degree  of  favoring  of  vibrational  energy  would  be  more  realistic.  Several  models  have  been 
implemented  and  analyzed  in  this  report.  All  models  allow  dissociation  (with  probability 
Pd)  only  of  a  diatom  with  the  collision  energy  ec  =  et  +  er  +  >  D.  (Pd  =  0  for  a 

diatom  with  insufficient  collision  energy  for  this  endothermic  reaction).  All  models  need 
to  be  calibrated  to  some  known  equilibrium  rate  constant  at  some  temperature  (or  some 
range  of  temperatures).  As  with  any  other  calibration  to  some  macroscopic  property,  the 
calibration  constant  will  be  dependent  on  the  details  of  the  collision  procedure.  Relatively 
comprehensive  comparisons  of  these  dissociation  models  have  been  made  and  are  documented 
as  part  of  this  report. 


1.  The  Exact  Available  Energy  (EAE)  model  of  Bird1  (Chap  6.7,  11.5,  see  also  Carlson 
&  Bird21),  where  dissociation  is  modeled  using  the  Borgnakke-Larsen  redistribution 
of  vibrational  energy,  and  a  diatom  dissociates  when  the  accepted  post-collision  vi¬ 
brational  energy  v'  >  vd-  The  model  requires  no  input  parameters,  but  a  calibration 
constant  is  allowed  for. 

2.  A  form  of  the  Weak  Vibrational  Bias  (WVB)  model,  given  by  Equation  (22)  of  Koura^ 


Pd  =  A; 


D  -  €y 


exp  (A[^  -  1]) 


(25) 


where  Ai  is  a  calibration  constant  chosen  to  allow  adjustment  of  the  overall  rate  con¬ 
stant,  and  A  is  an  input  parameter  controlling  the  degree  of  vibrational  favoring  or 
weighting.  The  symmetrized  cross  section  present  in  the  original  form  has  been  omit¬ 
ted  due  to  its  small  effect. 

3.  The  Full  Threshold  (FT)  model  of  Macheret,  et  al.22  The  model  characteristics  and 
appropriate  probabilities,  Equations  (6),  (10)  and  (16)-(19)  of  that  reference,  are  well 
described  in  that  report  and  are  not  repeated  here.  The  model  contains  no  adjustable 
parameters,  but  a  calibration  constant  Ai  has  been  added  here. 

4.  The  Simplified  Threshold  (ST)  model  of  Boyd,23  based  on  the  FT  model,  but  using 
simplified  forms  of  the  probability  equations  (Equations  (25)  and  (30)  in  Macheret,  et 
al.,  Equations  (4)- (9)  in  Boyd),  valid  only  for  the  limiting  case  of  high  translational 
energy  collisions,  et  >  D. 
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Exchange.  The  choice  of  the  model  for  exchange  is  dependent  on  the  assumed  degree  of 
vibrational  favoring  present  in  the  reaction.  As  shown  in  this  report,  the  dominant  exchange 
reactions  in  air,  the  Zeldovich  reactions,  are  not  (or  are  only  weakly)  vibrationally  favored. 
As  such,  the  EAE  model  as  described  by  Bird1  (Chap  6.8)  has  been  found  to  be  suitable 
for  describing  the  forward  (endothermic)  reaction  (in  the  sense  of  recovering  the  vibrational 
distribution  of  diatoms  chosen  to  undergo  exchange,  and  in  recovering  the  approximate 
Arrhenius  equilibrium  rate  constant).  This  model  is  implemented  analogously  to  that  for 
dissociation,  with  a  diatom  possessing  ec  =  et+er+ev  >  Ea,  (Ea  being  the  activation  energy) 
undergoing  exchange  when  the  sampled  post-collision  vibration  level  v'  >  vEa. 

Relatively  little  analysis  of  the  combined  forward  /  reverse  reaction  chemistry  model  has 
been  done  here.  Following  the  logic  of  Bird,  the  (potentially  temperature-dependent)  reverse 
reaction  probability  can  be  evaluated  using  the  equilibrium  rate  constant  and  the  calibrated 
form  of  the  forward  analytic  EAE  rate  constant,  evaluated  at  the  cell  temperature. 

Recombination.  For  the  low  gas  densities  of  interest,  recombination  reactions  (depen¬ 
dent  on  three-body  collisions)  are  negligible.  As  such,  relatively  little  analysis  of  recom- 
bination  models  has  been  done  here.  Models  of  varying  degrees  of  physical  validity  can 

be  generated  using  the  principle  of  detailed  balance  and  the  dissociation  models  described 
above. 

To  calibrate  the  combined  dissociation  /  recombination  models  to  some  known  equilib¬ 
rium  rate  (or  degree  of  dissociation),  limited  calculations  have  been  done  using  the  simple 
temperature-dependent  recombination  model  of  Bird1  (Chap  6.7,  11.5),  based  on  the  ana¬ 
lytic  EAE  dissociation  rate  constant.  Since  the  effective  rate  constant  extracted  from  the 
simulations  was  found  to  be  different  from  that  predicted  analytically,  the  recombination 
probability  must  be  recalibrated  to  recover  expected  the  equilibrium  constant  over  some 
range  of  temperature. 

Post-Collision  Redistribution  of  Energy.  A  particle  that  passes  the  acceptance  /  re¬ 
jection  test  (i.e.,  R  <  Px,  where  R  is  a  uniform  random  number)  shown  in  Figure  A2 
undergoes  redistribution  of  energy  in  that  mode  x  (for  the  case  of  internal  relaxation)  or 
undergoes  reaction  (for  the  case  of  chemical  reaction).  For  the  case  of  internal  relaxation, 
for  most  of  the  collision  models  considered  here,  the  redistribution  is  done  via  the  generalized 
Borgnakke-Larsen  procedure  (Bird,1  Chap.  5.5),  where  the  post-collision  value  is  sampled 
from  the  Borgnakke-Larsen  (equilibrium)  distribution  at  the  collision  energy  ec  —  et  +  ex, 
where  Q  denotes  the  relative  translational  energy  and  ex  denotes  the  energy  in  mode  x  of 
the  particle  undergoing  redistribution. 

Rotation.  Given  a  rotationally  inelastic  collision,  the  post-collision  rotational  level  j' 
is  sampled  using  a  modified  version  of  the  Boyd2  Borgnakke-Larsen  method.  The  peak- 
normalized  distribution  is  given  by  Equation  (10)  of  Boyd  as 

Q'(i'-e  )  =  W  +  Wtc-fU'  +  VMr]1-" 

’  (2jp  +  l)[ec  —  jp(jp  +  l)k9r]l~w 


(26) 


where  the  notation  of  Boyd  is  retained  for  u>(=  uB  -  1/2),  and  the  peak  occurs  at 


jp=  -1  + 


1  +  4  ec/kdr 


3 -2a ; 


l/2> 


/2 


(27) 


This  equation  differs  from  that  given  by  Equation  (11)  of  Boyd.  It  is,  however,  analogous 
to  that  given  by  Equation  (30)  of  Koura®  for  the  case  of  pair,  rather  than  single  particle, 
relaxation.  The  maximum  allowable  post-collision  jmax  is  limited  to  jo  if  ec  >  D;  otherwise 
it  is  as  given  by  Equation  (12)  of  Boyd  (the  present  Equation  (1),  with  D  replaced  by  ec). 
Comparing  jmax  and  jp,  since  ec  >  er,  the  relation  jmax  >  jp  will  always  be  satisfied,  and  the 
requirement  of  item  (3)  of  Boyd  is  irrelevant.  Note  that  in  this  method,  the  post-collision 
value  ?  is  allowed  to  take  any  value  0  <  j'  <  jmax,  including  j'  =  j,  and  also  that  there  is 
no  restriction  that  Sj  =  \j  —  j'\  be  even- valued  for  homonuclear  molecules.  These  constraints 
could  be  implemented  numerically,  however.  More  analysis  of  rotational  relaxation  will  be 
carried  out  in  the  future. 

Vibration.  The  present  analysis  includes  only  V-T  processes  and  ignores  the  physically- 
allowed  vibrational- vibrational  (V-V)  relaxation  process  in  diatom-diatom  collisions.  For  the 
Abe  state-specific  single-quantum  model,  the  post-collision  vibrational  level  is  determined 
directly  by  acceptance  (or  rejection)  of  Pv+i  =  PVtV+i/Pv  (or  Pv-\  —  1  —  Pt,+i).  A  similar 
method  is  used  for  the  Koura  ITFITS  formulation  (though  a  single  state  must  be  probabilis¬ 
tically  chosen  from  the  multiple-quanta  available). 

For  the  other  models,  v'  is  chosen  by  sampling  from  the  Borgnakke-Larsen  distribution 
as  discussed  by  Bird1  (Chap.  5.6),  where  0  <  v'  <  vmax,  with  ev(vmax)  <  ec.  Note  that  these 
models  allow  v'  =  v. 

Dissociation.  When  a  dissociation  reaction  occurs,  the  dissociation  energy  is  decre¬ 
mented  from  the  collision  energy  ec  =  et  +  er  +  e„  of  the  dissociating  diatom,  and  the 
remaining  energy  is  redistributed  to  the  relative  translational  energy  of  the  colliding  pair, 
and  then  to  the  relative  translational  energy  of  the  dissociating  atoms  following  the  method 
of  Bird1  (Chap  11.5). 

Exchange.  When  an  exchange  reaction  occurs,  the  heat  of  reaction  is  decremented  from 
the  collision  energy  ec  —  et  +  er  +  ev  of  the  colliding  diatom,  and  the  remaining  energy  is 
serially  redistributed  to  the  vibrational  and  rotational  mode  of  the  newly-created  diatom, 
and  then  to  the  relative  translational  energy  of  the  new  diatom-atom  pair. 

Recombination.  Following  the  Bird1  (Chap  11.5)  EAE  model,  the  new  diatom  is 
created  with  velocity  equal  to  the  center  of  mass  velocity  of  the  recombining  pair,  a  random 
third  body  is  chosen  for  collision  with  the  new  diatom,  and  the  energy  et  +  D  is  redistributed 
serially  into  the  vibrational  and  rotational  modes  of  the  newly-created  diatom,  and  then  to 
the  relative  translational  energy  of  the  diatom-third  body  pair. 
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Figure  A1 

Schematic  of  Serial  Collision  Procedure. 
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Figure  A2 

Schematic  of  Inelastic  or  Chemical  Reaction  Collision  Event  Selection 
Procedure  for  Diatom-Diatom  Pair,  with  Diatom  A2  chosen  first  (B2 

analogous). 
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